Commit e7344bf3 authored by Bruno Chareyre's avatar Bruno Chareyre
Browse files

- fix a bug that would let the capillary tables empty when postLoad was not triggered

parent 10022b85
......@@ -45,7 +45,7 @@ if model_type == 1:#hertz model with capillary forces
[Ip2_FrictMat_FrictMat_MindlinCapillaryPhys(label='ContactModel')],#for hertz model only
[Law2_ScGeom_MindlinPhys_Mindlin()]#for hertz model only
),
Law2_ScGeom_CapillaryPhys_Capillarity(CapillaryPressure=10000),#for hertz model only
Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=10000),#for hertz model only
NewtonIntegrator(damping=local_damping,gravity=(0,0,-9.81)),
]
ContactModel.betan=viscous_normal
......@@ -60,7 +60,7 @@ else:
[Ip2_FrictMat_FrictMat_CapillaryPhys()], #for linear model only
[Law2_ScGeom_FrictPhys_CundallStrack()], #for linear model only
),
Law2_ScGeom_CapillaryPhys_Capillarity(CapillaryPressure=10000),#for linear model only
Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=10000),#for linear model only
NewtonIntegrator(damping=local_damping,gravity=(0,0,-9.81)),
]
......
......@@ -15,16 +15,21 @@ class CapillaryPhys : public FrictPhys
virtual ~CapillaryPhys();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(CapillaryPhys,FrictPhys,"Physics (of interaction) for Law2_ScGeom_CapillaryPhys_Capillarity.",
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(CapillaryPhys,FrictPhys,"Physics (of interaction) for Law2_ScGeom_CapillaryPhys_Capillarity.",
((bool,meniscus,false,,"Presence of a meniscus if true"))
((bool,isBroken,false,,"If true, capillary force is zero and liquid bridge is inactive."))
((Real,CapillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
((Real,Vmeniscus,0.,,"Volume of the menicus"))
((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
((Real,vMeniscus,0.,,"Volume of the menicus"))
((Real,Delta1,0.,,"Defines the surface area wetted by the meniscus on the smallest grains of radius R1 (R1<R2)"))
((Real,Delta2,0.,,"Defines the surface area wetted by the meniscus on the biggest grains of radius R2 (R1<R2)"))
((Vector3r,Fcap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
((Vector3r,fCap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
((short int,fusionNumber,0.,,"Indicates the number of meniscii that overlap with this one"))
,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
,/*deprec*/
((Fcap,fCap,"naming convention"))
((CapillaryPressure,capillaryPressure,"naming convention"))
,,
createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
,
);
REGISTER_CLASS_INDEX(CapillaryPhys,FrictPhys);
};
......
......@@ -46,7 +46,7 @@ void CapillaryStressRecorder::action()
Real f1_cap_x=0, f1_cap_y=0, f1_cap_z=0, x1=0, y1=0, z1=0, x2=0, y2=0, z2=0;
Real sig11_cap=0, sig22_cap=0, sig33_cap=0, sig12_cap=0, sig13_cap=0,
sig23_cap=0, Vwater = 0, CapillaryPressure = 0;
sig23_cap=0, Vwater = 0, capillaryPressure = 0;
int j = 0;
InteractionContainer::iterator ii = scene->interactions->begin();
......@@ -67,7 +67,7 @@ void CapillaryStressRecorder::action()
unsigned int id1 = interaction -> getId1();
unsigned int id2 = interaction -> getId2();
Vector3r fcap = meniscusParameters->Fcap;
Vector3r fcap = meniscusParameters->fCap;
f1_cap_x=fcap[0];
f1_cap_y=fcap[1];
......@@ -94,8 +94,8 @@ void CapillaryStressRecorder::action()
/// Liquid volume
Vwater += meniscusParameters->Vmeniscus;
CapillaryPressure=meniscusParameters->CapillaryPressure;
Vwater += meniscusParameters->vMeniscus;
capillaryPressure=meniscusParameters->capillaryPressure;
}
......@@ -148,7 +148,7 @@ void CapillaryStressRecorder::action()
<< lexical_cast<string>(SIG_12_cap) << " "
<< lexical_cast<string>(SIG_13_cap)<< " "
<< lexical_cast<string>(SIG_23_cap)<< " "
<< lexical_cast<string>(CapillaryPressure) << " "
<< lexical_cast<string>(capillaryPressure) << " "
<< lexical_cast<string>(Sr)<< " "
<< lexical_cast<string>(w)<< " "
<< endl;
......
......@@ -338,7 +338,7 @@ void CapillaryTriaxialTest::createActors(shared_ptr<Scene>& scene)
// capillary
shared_ptr<Law2_ScGeom_CapillaryPhys_Capillarity> capillaryCohesiveLaw(new Law2_ScGeom_CapillaryPhys_Capillarity);
capillaryCohesiveLaw->CapillaryPressure = CapillaryPressure;
capillaryCohesiveLaw->capillaryPressure = capillaryPressure;
// capillaryCohesiveLaw->fusionDetection = fusionDetection;
// capillaryCohesiveLaw->binaryFusion = binaryFusion;
......
......@@ -71,7 +71,7 @@ class CapillaryTriaxialTest : public FileGenerator
((string,importFilename,"",,"File with positions and sizes of spheres."))
((string,Key,"",,"A code that is added to output filenames."))
((string,fixedBoxDims,"",,"string that contains some subset (max. 2) of {'x','y','z'} ; contains axes will have box dimension hardcoded, even if box is scaled as mean_radius is prescribed: scaling will be applied on the rest."))
((Real,CapillaryPressure,0,,"Define succion in the packing [Pa]. This is the value used in the capillary model."))
((Real,capillaryPressure,0,,"Define succion in the packing [Pa]. This is the value used in the capillary model."))
((bool,water,true,,"activate capillary model"))
((bool,fusionDetection,false,,"test overlaps between liquid bridges on modify forces if overlaps exist"))
((bool,binaryFusion,true,,"Defines how overlapping bridges affect the capillary forces (see :yref:`CapillaryTriaxialTest::fusionDetection`). If binary=true, the force is null as soon as there is an overlap detected, if not, the force is divided by the number of overlaps."))
......
......@@ -19,16 +19,20 @@ class MindlinCapillaryPhys : public MindlinPhys
virtual ~MindlinCapillaryPhys();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(MindlinCapillaryPhys,MindlinPhys,"Adds capillary physics to Mindlin's interaction physics.",
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(MindlinCapillaryPhys,MindlinPhys,"Adds capillary physics to Mindlin's interaction physics.",
((bool,meniscus,false,,"Presence of a meniscus if true"))
((bool,isBroken,false,,"If true, capillary force is zero and liquid bridge is inactive."))
((Real,CapillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
((Real,Vmeniscus,0.,,"Volume of the menicus"))
((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
((Real,vMeniscus,0.,,"Volume of the menicus"))
((Real,Delta1,0.,,"Defines the surface area wetted by the meniscus on the smallest grains of radius R1 (R1<R2)"))
((Real,Delta2,0.,,"Defines the surface area wetted by the meniscus on the biggest grains of radius R2 (R1<R2)"))
((Vector3r,Fcap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
((Vector3r,fCap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
((short int,fusionNumber,0.,,"Indicates the number of meniscii that overlap with this one"))
,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
,/*deprec*/
((Fcap,fCap,"naming convention"))
((CapillaryPressure,capillaryPressure,"naming convention"))
,,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
,
);
REGISTER_CLASS_INDEX(MindlinCapillaryPhys,MindlinPhys);
};
......
......@@ -67,6 +67,7 @@ MeniscusParameters::~MeniscusParameters()
void Law2_ScGeom_CapillaryPhys_Capillarity::action()
{
if (!scene) cerr << "scene not defined!";
if (!capillary) postLoad(*this);//when the script does not define arguments, postLoad is never called
shared_ptr<BodyContainer>& bodies = scene->bodies;
if (fusionDetection && !bodiesMenisciiList.initialized) bodiesMenisciiList.prepare(scene);
......@@ -119,7 +120,7 @@ void Law2_ScGeom_CapillaryPhys_Capillarity::action()
Real D = alpha*((b2->state->pos-b1->state->pos).norm()-(currentContactGeometry->radius1+ currentContactGeometry->radius2)); // scGeom->penetrationDepth could probably be used here?
if ((currentContactGeometry->penetrationDepth>=0)|| D<=0 || createDistantMeniscii) { //||(scene->iter < 1) ) // a simplified way to define meniscii everywhere
D=0; // defines Fcap when spheres interpenetrate. D<0 leads to wrong interpolation has D<0 has no solution in the interpolation : this is not physically interpretable!! even if, interpenetration << grain radius.
D=0; // defines fCap when spheres interpenetrate. D<0 leads to wrong interpolation has D<0 has no solution in the interpolation : this is not physically interpretable!! even if, interpenetration << grain radius.
if (!hertzOn) {
if (fusionDetection && !cundallContactPhysics->meniscus) bodiesMenisciiList.insert((*ii));
cundallContactPhysics->meniscus=true;
......@@ -132,10 +133,10 @@ void Law2_ScGeom_CapillaryPhys_Capillarity::action()
/// Suction (Capillary pressure):
Real Pinterpol = 0;
if (!hertzOn) Pinterpol = cundallContactPhysics->isBroken ? 0 : CapillaryPressure*(R2/liquidTension);
else Pinterpol = mindlinContactPhysics->isBroken ? 0 : CapillaryPressure*(R2/liquidTension);
if (!hertzOn) cundallContactPhysics->CapillaryPressure = CapillaryPressure;
else mindlinContactPhysics->CapillaryPressure = CapillaryPressure;
if (!hertzOn) Pinterpol = cundallContactPhysics->isBroken ? 0 : capillaryPressure*(R2/liquidTension);
else Pinterpol = mindlinContactPhysics->isBroken ? 0 : capillaryPressure*(R2/liquidTension);
if (!hertzOn) cundallContactPhysics->capillaryPressure = capillaryPressure;
else mindlinContactPhysics->capillaryPressure = capillaryPressure;
/// Capillary solution finder:
if ((Pinterpol>=0) && (hertzOn? mindlinContactPhysics->meniscus : cundallContactPhysics->meniscus)) {
......@@ -145,17 +146,17 @@ void Law2_ScGeom_CapillaryPhys_Capillarity::action()
solution(Pinterpol? capillary->Interpolate(R1,R2,Dinterpol, Pinterpol, currentIndexes) : MeniscusParameters());
/// capillary adhesion force
Real Finterpol = solution.F;
Vector3r Fcap = Finterpol*(2*Mathr::PI*(R2/alpha)*liquidTension)*currentContactGeometry->normal;
if (!hertzOn) cundallContactPhysics->Fcap = Fcap;
else mindlinContactPhysics->Fcap = Fcap;
Vector3r fCap = Finterpol*(2*Mathr::PI*(R2/alpha)*liquidTension)*currentContactGeometry->normal;
if (!hertzOn) cundallContactPhysics->fCap = fCap;
else mindlinContactPhysics->fCap = fCap;
/// meniscus volume
Real Vinterpol = solution.V * pow(R2/alpha,3);
if (!hertzOn) {
cundallContactPhysics->Vmeniscus = Vinterpol;
cundallContactPhysics->vMeniscus = Vinterpol;
if (Vinterpol != 0) cundallContactPhysics->meniscus = true;
else cundallContactPhysics->meniscus = false;
} else {
mindlinContactPhysics->Vmeniscus = Vinterpol;
mindlinContactPhysics->vMeniscus = Vinterpol;
if (Vinterpol != 0) mindlinContactPhysics->meniscus = true;
else mindlinContactPhysics->meniscus = false;
}
......@@ -190,15 +191,15 @@ void Law2_ScGeom_CapillaryPhys_Capillarity::action()
short int& fusionNumber = hertzOn?mindlinContactPhysics->fusionNumber:cundallContactPhysics->fusionNumber;
if (binaryFusion) {
if (fusionNumber!=0) { //cerr << "fusion" << endl;
hertzOn?mindlinContactPhysics->Fcap:cundallContactPhysics->Fcap = Vector3r::Zero();
hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap = Vector3r::Zero();
continue;
}
}
//LINEAR VERSION : capillary force is divided by (fusionNumber + 1) - NOTE : any decreasing function of fusionNumber can be considered in fact
else if (fusionNumber !=0) hertzOn?mindlinContactPhysics->Fcap:cundallContactPhysics->Fcap /= (fusionNumber+1.);
else if (fusionNumber !=0) hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap /= (fusionNumber+1.);
}
scene->forces.addForce((*ii)->getId1(), hertzOn?mindlinContactPhysics->Fcap:cundallContactPhysics->Fcap);
scene->forces.addForce((*ii)->getId2(),-(hertzOn?mindlinContactPhysics->Fcap:cundallContactPhysics->Fcap));
scene->forces.addForce((*ii)->getId1(), hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap);
scene->forces.addForce((*ii)->getId2(),-(hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap));
}
}
}
......@@ -358,7 +359,7 @@ Tableau::Tableau(const char* filename)
for (int i=0; i<n_D; i++)
full_data.push_back(TableauD(file));
file.close();
//cerr << *this; // exemple d'utilisation de la fonction d'ecriture (this est le pointeur vers l'objet courant)
}
Tableau::~Tableau()
......
......@@ -13,7 +13,6 @@
#include <yade/core/GlobalEngine.hpp>
#include <set>
#include <boost/tuple/tuple.hpp>
#include <vector>
#include <list>
#include <utility>
......@@ -94,12 +93,15 @@ class Law2_ScGeom_CapillaryPhys_Capillarity : public GlobalEngine
void action();
void postLoad(Law2_ScGeom_CapillaryPhys_Capillarity&);
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CapillaryPhys_Capillarity,GlobalEngine,"This law allows one to take into account capillary forces/effects between spheres coming from the presence of interparticular liquid bridges (menisci).\n\nThe control parameter is the capillary pressure (or suction) Uc = ugas - Uliquid. Liquid bridges properties (volume V, extent over interacting grains delta1 and delta2) are computed as a result of the defined capillary pressure and of the interacting geometry (spheres radii and interparticular distance).\n\nReferences: in english [Scholtes2009b]_; more detailed, but in french [Scholtes2009d]_.\n\nThe law needs ascii files M(r=i) with i=R1/R2 to work (see https://yade-dem.org/index.php/CapillaryTriaxialTest). These ASCII files contain a set of results from the resolution of the Laplace-Young equation for different configurations of the interacting geometry.\n\nIn order to allow capillary forces between distant spheres, it is necessary to enlarge the bounding boxes using :yref:`Bo1_Sphere_Aabb::aabbEnlargeFactor` and make the Ig2 define define distant interactions via :yref:`interactionDetectionFactor<Ig2_Sphere_Sphere_ScGeom::interactionDetectionFactor>`. It is also necessary to disable interactions removal by the constitutive law (:yref:`Law2<Law2_ScGeom_FrictPhys_CundallStrack::neverErase>=True`). The only combinations of laws supported are currently capillary law + :yref:`Law2_ScGeom_FrictPhys_CundallStrack` and capillary law + :yref:`Law2_ScGeom_MindlinPhys_Mindlin` (and the other variants of Hertz-Mindlin).\n\nSee CapillaryPhys-example.py for an example script.",
((Real,CapillaryPressure,0.,,"Value of the capillary pressure Uc defines as Uc=Ugas-Uliquid"))
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(Law2_ScGeom_CapillaryPhys_Capillarity,GlobalEngine,"This law allows one to take into account capillary forces/effects between spheres coming from the presence of interparticular liquid bridges (menisci).\n\nThe control parameter is the capillary pressure (or suction) Uc = ugas - Uliquid. Liquid bridges properties (volume V, extent over interacting grains delta1 and delta2) are computed as a result of the defined capillary pressure and of the interacting geometry (spheres radii and interparticular distance).\n\nReferences: in english [Scholtes2009b]_; more detailed, but in french [Scholtes2009d]_.\n\nThe law needs ascii files M(r=i) with i=R1/R2 to work (see https://yade-dem.org/index.php/CapillaryTriaxialTest). These ASCII files contain a set of results from the resolution of the Laplace-Young equation for different configurations of the interacting geometry.\n\nIn order to allow capillary forces between distant spheres, it is necessary to enlarge the bounding boxes using :yref:`Bo1_Sphere_Aabb::aabbEnlargeFactor` and make the Ig2 define define distant interactions via :yref:`interactionDetectionFactor<Ig2_Sphere_Sphere_ScGeom::interactionDetectionFactor>`. It is also necessary to disable interactions removal by the constitutive law (:yref:`Law2<Law2_ScGeom_FrictPhys_CundallStrack::neverErase>=True`). The only combinations of laws supported are currently capillary law + :yref:`Law2_ScGeom_FrictPhys_CundallStrack` and capillary law + :yref:`Law2_ScGeom_MindlinPhys_Mindlin` (and the other variants of Hertz-Mindlin).\n\nSee CapillaryPhys-example.py for an example script.",
((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Uc=Ugas-Uliquid"))
((bool,fusionDetection,false,,"If true potential menisci overlaps are checked"))
((bool,binaryFusion,true,,"If true, capillary forces are set to zero as soon as, at least, 1 overlap (menisci fusion) is detected"))
((bool,hertzOn,false,,"|yupdate| true if hertz model is used"))
((bool,createDistantMeniscii,false,,"Generate meniscii between distant spheres? Else only maintain the existing one. For modeling a wetting path this flag should always be false. For a drying path it should be true for one step (initialization) then false, as in the logic of [Scholtes2009c]_"))
((bool,createDistantMeniscii,false,,"Generate meniscii between distant spheres? Else only maintain the existing ones. For modeling a wetting path this flag should always be false. For a drying path it should be true for one step (initialization) then false, as in the logic of [Scholtes2009c]_"))
,/*deprec*/
((CapillaryPressure,capillaryPressure,"naming convention"))
,,,
);
};
......
......@@ -60,12 +60,12 @@ void SampleCapillaryPressureEngine::action()
if (scene->iter % 100 == 0) cerr << "pressure variation!!" << endl;
if ((Pressure>=0) && (Pressure<=1000000000)) Pressure += PressureVariation;
capillaryCohesiveLaw->CapillaryPressure = Pressure;
capillaryCohesiveLaw->capillaryPressure = Pressure;
capillaryCohesiveLaw->fusionDetection = fusionDetection;
capillaryCohesiveLaw->binaryFusion = binaryFusion;
}
else { capillaryCohesiveLaw->CapillaryPressure = Pressure;
else { capillaryCohesiveLaw->capillaryPressure = Pressure;
capillaryCohesiveLaw->fusionDetection = fusionDetection;
capillaryCohesiveLaw->binaryFusion = binaryFusion;}
if (scene->iter % 100 == 0) cerr << "capillary pressure = " << Pressure << endl;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment