------------------------------------------------------------ revno: 3707 committer: Chao Yuan <chaoyuan2...@gmail.com> timestamp: Thu 2015-07-23 19:26:21 +0200 message: add a new version of capillary law. (pushed by Caroline) added: pkg/dem/CapillaryPhys1.cpp pkg/dem/CapillaryPhys1.hpp pkg/dem/DelaunayInterpolation.hpp pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.cpp pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.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
=== added file 'pkg/dem/CapillaryPhys1.cpp' --- pkg/dem/CapillaryPhys1.cpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/CapillaryPhys1.cpp 2015-07-23 17:26:21 +0000 @@ -0,0 +1,58 @@ +//keep this #ifdef as long as you don't really want to realize a final version publicly, it will save compilation time for everyone else +//when you want it compiled, you can just uncomment the following line +// #define CAPILLARYPHYS1 +#ifdef CAPILLARYPHYS1 + +#include <pkg/dem/CapillaryPhys1.hpp> +#include<pkg/dem/ScGeom.hpp> + +#include<core/Omega.hpp> +#include<core/Scene.hpp> + +CapillaryPhys1::~CapillaryPhys1() +{ +} + +YADE_PLUGIN((CapillaryPhys1)); + + + +YADE_PLUGIN((Ip2_FrictMat_FrictMat_CapillaryPhys1)); + +void Ip2_FrictMat_FrictMat_CapillaryPhys1::go( const shared_ptr<Material>& b1 //FrictMat + , const shared_ptr<Material>& b2 // FrictMat + , const shared_ptr<Interaction>& interaction) +{ + ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get()); + if(geom) + { + + if(!interaction->phys) + { + const shared_ptr<FrictMat>& sdec1 = YADE_PTR_CAST<FrictMat>(b1); + const shared_ptr<FrictMat>& sdec2 = YADE_PTR_CAST<FrictMat>(b2); + + if (!interaction->phys) interaction->phys = shared_ptr<CapillaryPhys1>(new CapillaryPhys1()); + const shared_ptr<CapillaryPhys1>& contactPhysics = YADE_PTR_CAST<CapillaryPhys1>(interaction->phys); + + Real Ea = sdec1->young; + Real Eb = sdec2->young; + Real Va = sdec1->poisson; + Real Vb = sdec2->poisson; + Real Da = geom->radius1; // FIXME - multiply by factor of sphere interaction distance (so sphere interacts at bigger range that its geometrical size) + Real Db = geom->radius2; // FIXME - as above + Real fa = sdec1->frictionAngle; + Real fb = sdec2->frictionAngle; + Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses + Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere + + contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb)); + contactPhysics->kn = Kn; + contactPhysics->ks = Ks; + contactPhysics->computeBridge =computeDefault; + } + } +}; + +#endif //CAPILLARYPHYS1 + === added file 'pkg/dem/CapillaryPhys1.hpp' --- pkg/dem/CapillaryPhys1.hpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/CapillaryPhys1.hpp 2015-07-23 17:26:21 +0000 @@ -0,0 +1,99 @@ +/************************************************************************* +* Copyright (C) 2006 by luc Scholtes * +* luc.schol...@hmg.inpg.fr * +* * +* This program is free software; it is licensed under the terms of the * +* GNU General Public License v2 or later. See file LICENSE for details. * +*************************************************************************/ +#pragma once +#include<pkg/dem/FrictPhys.hpp> +#include<pkg/dem/DelaunayInterpolation.hpp> +#include<pkg/common/Dispatching.hpp> +#include<pkg/common/ElastMat.hpp> + +class MeniscusPhysicalData { +public: + double R; + double volume; + double distance; + double surface; + double energy; + double force; + double succion; + double delta1; + double delta2; + double arcLength; + bool ending; + //default ctor + MeniscusPhysicalData() : R(0), volume(0), distance(0), surface(0), energy(0), force(0), succion(0), delta1(0), delta2(0), arcLength(0), ending(1) {} + //ctor with input values + MeniscusPhysicalData(const double& r, const double& v, const double& d, const double& s, const double& e, const double& f, const double& p, const double& a1, const double& a2, const double& arc, bool end) : R(r), volume(v), distance(d), surface(s), energy(e), force(f), succion(p), delta1(a1), delta2(a2), arcLength(arc), ending(end) {} + + //a minimal list of operators for the interpolation + //these operators are requirements for the DataType template parameter of interpolate() + //FIXME: algebra operations must include energy, perimeter, and any other new variable for including them in the interpolation + MeniscusPhysicalData& operator+= (const MeniscusPhysicalData& m2) { + R+=m2.R; volume+=m2.volume; distance+=m2.distance; surface+=m2.surface; energy+=m2.energy; force+=m2.force; succion+=m2.succion; delta1+=m2.delta1; delta2+=m2.delta2; arcLength+=m2.arcLength; ending=(ending && m2.ending); + return *this;} + + MeniscusPhysicalData operator* (const double& fact) const { + return MeniscusPhysicalData(fact*R, fact*volume, fact*distance, fact*surface, fact*energy, fact*force, fact*succion, fact*delta1, fact*delta2, fact*arcLength, ending);} + + const MeniscusPhysicalData& operator= (const MeniscusPhysicalData& m1) { + R=m1.R; volume=m1.volume; distance=m1.distance; surface=m1.surface; energy=m1.energy; force=m1.force; succion=m1.succion; delta1=m1.delta1; delta2=m1.delta2; arcLength=m1.arcLength; ending=m1.ending; + return *this;} +}; + +//The structure for the meniscus: physical properties + cached values for fast interpolation +class Meniscus { + public: + typedef MeniscusPhysicalData Data; + Data data; //the variables of Laplace's problem + DT::Cell_handle cell; //pointer to the last location in the triangulation, for faster locate() + std::vector<K::Vector_3> normals;// 4 normals relative to the current cell + + Meniscus() : data(), cell(DT::Cell_handle()), normals(std::vector<K::Vector_3>()) {} +}; + +class CapillaryPhys1 : public FrictPhys +{ + public : + int currentIndexes [4]; // used for faster interpolation (stores previous positions in tables) + Meniscus m; + + virtual ~CapillaryPhys1(); + + YADE_CLASS_BASE_DOC_ATTRS_CTOR(CapillaryPhys1,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.")) + ((bool,computeBridge,true,,"If true, capillary bridge will be computed if not it will be ignored.")) + ((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")) + ((Real,SInterface,0.,,"Fluid-Gaz Interfacial area")) + ((Real,arcLength,0.,,"Arc Length of the Fluid-Gaz Interface")) + ((short int,fusionNumber,0.,,"Indicates the number of meniscii that overlap with this one")) + ,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0; + ); + REGISTER_CLASS_INDEX(CapillaryPhys1,FrictPhys); +}; +REGISTER_SERIALIZABLE(CapillaryPhys1); + +class Ip2_FrictMat_FrictMat_CapillaryPhys1 : public IPhysFunctor +{ + public : + virtual void go( const shared_ptr<Material>& b1, + const shared_ptr<Material>& b2, + const shared_ptr<Interaction>& interaction); + + FUNCTOR2D(FrictMat,FrictMat); + YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_FrictMat_FrictMat_CapillaryPhys1,IPhysFunctor, "RelationShips to use with Law2_ScGeom_CapillaryPhys_Capillarity1\n\n In these RelationShips all the interaction attributes are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.", + ((bool,computeDefault,true,,"bool to assign the default value of computeBridge.")),; + + ); + +}; +REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_CapillaryPhys1); + === added file 'pkg/dem/DelaunayInterpolation.hpp' --- pkg/dem/DelaunayInterpolation.hpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/DelaunayInterpolation.hpp 2015-07-23 17:26:21 +0000 @@ -0,0 +1,150 @@ +/************************************************************************* +* Copyright (C) 2013 by Bruno Chareyre <bruno.chare...@hmg.inpg.fr> * +* * +* This program is free software; it is licensed under the terms of the * +* GNU General Public License v2 or later. See file LICENSE for details. * +*************************************************************************/ + +#pragma once +#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> +#include <CGAL/Delaunay_triangulation_3.h> +#include <CGAL/Triangulation_vertex_base_with_info_3.h> +#include <CGAL/Cartesian.h> +#include <iostream> +#include <fstream> +///////////////////////////////////////////////////////////////////// +/* +TYPES: +Triangulation_vertex_base_with_id_3: we redefine a vertex base including an index for each vertex (available in CGAL in 2D but not in 3D), +MeniscusPhysicalData: the physical variables describing a capillary bridge, with a few algebraic operators for computing weighted averages +Meniscus: a structure combining MeniscusPhysicalData with some cached data allowing faster operations in multiple queries (pointer to the last cell found and its normals) + +FUNCTIONS: +getIncidentVtxWeights: an interpolation algorithm which is 20x faster than the natural neighbor interpolation of CGAL. + returns a list of vertices incident to a query point and their respective weights +interpolate: uses the results of getIncidentVtxWeights combined with a data array to return a weighted average, + may be used with arbitrary data types provided they have the required algebraic operators +main: example usage +*/ +///////////////////////////////////////////////////////////////////// + +namespace CGAL { + + + +//Vertex base including an index for each vertex, adapted from CGAL::Triangulation_vertex_base_with_id_2 +template < typename GT, typename Vb = Triangulation_vertex_base_with_info_3<unsigned,GT> > +class Triangulation_vertex_base_with_id_3 : public Vb +{ + int _id; +public: + typedef typename Vb::Cell_handle Cell_handle; + typedef typename Vb::Point Point; + template < typename TDS3 > + struct Rebind_TDS { + typedef typename Vb::template Rebind_TDS<TDS3>::Other Vb3; + typedef Triangulation_vertex_base_with_id_3<GT, Vb3> Other; + }; + + Triangulation_vertex_base_with_id_3(): Vb() {} + Triangulation_vertex_base_with_id_3(const Point & p): Vb(p) {} + Triangulation_vertex_base_with_id_3(const Point & p, Cell_handle c): Vb(p, c) {} + Triangulation_vertex_base_with_id_3(Cell_handle c): Vb(c) {} + unsigned int id() const { return this->info(); } + unsigned int& id() { return this->info(); } +}; + +// The function returning vertices and their weights for an arbitrary point in R3 space. +// The returned triplet contains: +// - the output iterator pointing to the filled vector of vertices +// - the sum of weights for normalisation +// - a bool telling if the query point was located inside the convex hull of the data points +template <class Dt, class OutputIterator> +Triple< OutputIterator, // iterator with value type std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT> + typename Dt::Geom_traits::FT, // Should provide 0 and 1 + bool > +getIncidentVtxWeights(const Dt& dt, + const typename Dt::Geom_traits::Point_3& Q, + OutputIterator nn_out, typename Dt::Geom_traits::FT & norm_coeff, + std::vector<typename Dt::Geom_traits::Vector_3>& normals, + typename Dt::Cell_handle& start = CGAL_TYPENAME_DEFAULT_ARG Dt::Cell_handle()) +{ + //helpful array for permutations + const int comb [6] = {1, 2, 3, 0, 1, 2}; + typedef typename Dt::Geom_traits Gt; + typedef typename Gt::Point_3 Point; + typedef typename Dt::Cell_handle Cell_handle; + typedef typename Dt::Locate_type Locate_type; + typedef typename Gt::FT Coord_type; + CGAL_triangulation_precondition (dt.dimension()== 3); + Locate_type lt; int li, lj; + Cell_handle c = dt.locate( Q, lt, li, lj, start); + bool updateNormals = (c!=start || normals.size()<4); + if (updateNormals) normals.clear(); + if ( lt == Dt::VERTEX ) + { + *nn_out++= std::make_pair(c->vertex(li),Coord_type(1)); + return make_triple(nn_out,norm_coeff=Coord_type(1),true); + } + else if (dt.is_infinite(c)) + return make_triple(nn_out, Coord_type(1), false);//point outside the convex-hull + norm_coeff=0; + for ( int k=0;k<4;k++ ) + { + if (updateNormals) { + normals.push_back(cross_product(c->vertex(comb[k])->point()-c->vertex(comb[k+1])->point(),c->vertex(comb[k])->point()-c->vertex(comb[k+2])->point())); + normals[k] = normals[k]/ + ((c->vertex(k)->point()-c->vertex(comb[k])->point())*normals[k]); + } + Coord_type closeness = ((Q-c->vertex(comb[k])->point())*normals[k]); + Coord_type w = closeness; + *nn_out++= std::make_pair(c->vertex(k),w); + norm_coeff += w; + } + start = c; + return make_triple(nn_out,norm_coeff,true); +} + + + +} //END NAMESPACE CGAL + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Delaunay_triangulation_3<K>::Geom_traits Traits; +typedef CGAL::Triangulation_vertex_base_with_id_3<Traits> Vb; +typedef CGAL::Triangulation_cell_base_3<Traits> Cb; +typedef CGAL::Triangulation_data_structure_3<Vb, Cb> Tds; +typedef CGAL::Delaunay_triangulation_3<Traits,Tds> DT; +typedef std::vector< std::pair< DT::Vertex_handle, K::FT> > Vertex_weight_vector; + +template <class Dt, class DataOwner> +typename DataOwner::Data interpolate1 (const Dt& dt, const typename Dt::Geom_traits::Point_3& Q, DataOwner& owner, const std::vector<typename DataOwner::Data>& rawData, bool reset) +{ + K::FT norm; + Vertex_weight_vector coords; + if (reset) owner.cell = dt.cells_begin(); + CGAL::Triple<std::back_insert_iterator<Vertex_weight_vector>,K::FT, bool> result = CGAL::getIncidentVtxWeights(dt, Q,std::back_inserter(coords), norm, owner.normals , owner.cell); + + typename DataOwner::Data data = typename DataOwner::Data();//initialize null solution + if (!result.third) return data;// out of the convex hull, we return the null solution + //else, we compute the weighted sum + for (unsigned int k=0; k<coords.size(); k++) data += (rawData[coords[k].first->id()]*coords[k].second); + if (!data.ending) return data*(1./result.second); + else return typename DataOwner::Data(); +} +template <class Dt, class DataOwner> +typename DataOwner::Data interpolate2 (const Dt& dt, const typename Dt::Geom_traits::Point_3& Q, DataOwner& owner, const std::vector<typename DataOwner::Data>& rawData, bool reset) +{ + K::FT norm; + Vertex_weight_vector coords; + if (reset) owner.cell = dt.cells_begin(); + CGAL::Triple<std::back_insert_iterator<Vertex_weight_vector>,K::FT, bool> result = CGAL::getIncidentVtxWeights(dt, Q,std::back_inserter(coords), norm, owner.normals , owner.cell); + + typename DataOwner::Data data = typename DataOwner::Data();//initialize null solution + if (!result.third) return data;// out of the convex hull, we return the null solution + //else, we compute the weighted sum + for (unsigned int k=0; k<coords.size(); k++) data += (rawData[coords[k].first->id()]*coords[k].second); + if (!data.ending) return data*(1./result.second); + else return typename DataOwner::Data(); +} + === added file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.cpp' --- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.cpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.cpp 2015-07-23 17:26:21 +0000 @@ -0,0 +1,606 @@ +/************************************************************************* +* Copyright (C) 2006 by luc Scholtes * +* luc.schol...@hmg.inpg.fr * +* * +* This program is free software; it is licensed under the terms of the * +* GNU General Public License v2 or later. See file LICENSE for details. * +*************************************************************************/ + +//Modifs : Parameters renamed as MeniscusParameters +//id1/id2 as id1 is the smallest grain, FIXME : wetting angle? +//FIXME : in triaxialStressController, change test about null force in updateStiffnessccc +//keep this #ifdef as long as you don't really want to realize a final version publicly, it will save compilation time for everyone else +//when you want it compiled, you can just uncomment the following line +// #define LAW2_SCGEOM_CAPILLARYPHYS_Capillarity1 +#ifdef LAW2_SCGEOM_CAPILLARYPHYS_Capillarity1 + +#include <pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.hpp> +#include <pkg/common/ElastMat.hpp> + +#include <pkg/dem/ScGeom.hpp> +#include <pkg/dem/HertzMindlin.hpp> +#include <core/Omega.hpp> +#include <core/Scene.hpp> +#include <lib/base/Math.hpp> +#include <iostream> +#include <fstream> + + +DT Law2_ScGeom_CapillaryPhys_Capillarity1::dtVbased; +DT Law2_ScGeom_CapillaryPhys_Capillarity1::dtPbased; + + Real Law2_ScGeom_CapillaryPhys_Capillarity1::intEnergy() +{ + Real energy=0; + FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){ + if(!I->isReal()) continue; + ScGeom* currentGeometry = static_cast<ScGeom*>(I->geom.get()); + CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>(I->phys.get()); + if(phys) { + if (phys->SInterface!=0){ + energy += liquidTension*(phys->SInterface-4*3.141592653589793238462643383279502884*(pow(currentGeometry->radius1,2))-4*3.141592653589793238462643383279502884*(pow(currentGeometry->radius2,2)));} + } + } + return energy; +} + + + Real Law2_ScGeom_CapillaryPhys_Capillarity1::wnInterface() +{ + Real wn=0; + FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){ + if(!I->isReal()) continue; + ScGeom* currentGeometry = static_cast<ScGeom*>(I->geom.get()); + CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>(I->phys.get()); + if(phys) { + if (phys->SInterface!=0){ + + wn += (phys->SInterface-2*3.141592653589793238462643383279502884*(pow(currentGeometry->radius1,2)*(1+cos(phys->Delta1))+pow(currentGeometry->radius2,2)*(1+cos(phys->Delta2)))); + } + } + } + return wn; +} + + Real Law2_ScGeom_CapillaryPhys_Capillarity1::swInterface() +{ + Real sw=0; + FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){ + if(!I->isReal()) continue; + ScGeom* currentGeometry = static_cast<ScGeom*>(I->geom.get()); + CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>(I->phys.get()); + if(phys) { + sw += (2*3.141592653589793238462643383279502884*(pow(currentGeometry->radius1,2)*(1-cos(phys->Delta1))+pow(currentGeometry->radius2,2)*(1-cos(phys->Delta2)))); + + } + } + return sw; +} + + + Real Law2_ScGeom_CapillaryPhys_Capillarity1::waterVolume() +{ + Real volume=0; + FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){ + if(!I->isReal()) continue; + CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>(I->phys.get()); + if(phys) { + volume += phys->vMeniscus;} + } + return volume; +} + +void Law2_ScGeom_CapillaryPhys_Capillarity1::triangulateData() { + /// We get data from a file and input them in triangulations + if (solutions.size()>0) {LOG_WARN("Law2_ScGeom_CapillaryPhys_Capillarity1 asking triangulation for the second time. Ignored."); return;} + ifstream file (inputFilename.c_str()); + if (!file.is_open()) { LOG_ERROR("No data file found for capillary law. Check path and inputFilename."); return;} + + // convention R,v,d,s,e,f,p,a1,a2,dummy (just for the example, define your own, + // dummy is because has too much values per line - with one useless extra colum,) + MeniscusPhysicalData dat; + double ending; + while ( file.good() ) { + file >>dat.succion>>dat.force>>dat.distance>>dat.volume>>dat.surface>>dat.arcLength>>dat.delta1>>dat.delta2>>dat.R>>ending; + dat.ending=(bool) ending; + solutions.push_back(dat); + } + file.close(); + // Make lists of points with index, so we can use range insertion, more efficient + // see http://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation_3SettingInformationWhileInserting + std::vector< std::pair<K::Point_3,unsigned> > pointsP, pointsV; + for (unsigned int k=0; k<solutions.size(); k++) { + pointsP.push_back(std::make_pair(K::Point_3(solutions[k].R, solutions[k].succion, solutions[k].distance),k)); + pointsV.push_back(std::make_pair(K::Point_3(solutions[k].R, solutions[k].volume, solutions[k].distance),k)); + } + // and now range insertion + dtPbased.insert(pointsP.begin(), pointsP.end()); + dtVbased.insert(pointsV.begin(), pointsV.end()); +} + +YADE_PLUGIN((Law2_ScGeom_CapillaryPhys_Capillarity1)); + +int firstIteration=1; +int x=0; + +void Law2_ScGeom_CapillaryPhys_Capillarity1::action() +{ + bool switched = (switchTriangulation == (imposePressure or totalVolumeConstant)); + switchTriangulation = (imposePressure or totalVolumeConstant); + InteractionContainer::iterator ii = scene->interactions->begin(); + InteractionContainer::iterator iiEnd = scene->interactions->end(); + if (imposePressure) { + solver(capillaryPressure,switched); + } + else{ + if (((totalVolumeConstant || (!totalVolumeConstant && firstIteration==1)) && totalVolumeofWater!=-1) || (totalVolumeConstant && totalVolumeofWater==-1)) + { + if (!totalVolumeConstant) x=1; + totalVolumeConstant=1; + Real p0=capillaryPressure; + Real slope; + Real eps=0.0000000001; + solver(p0,switched); + Real V0=waterVolume(); + if (totalVolumeConstant && totalVolumeofWater==-1 && firstIteration==1){ + totalVolumeofWater=V0; + firstIteration+=1; + } + Real p1=capillaryPressure+0.1; + solver(p1,switched); + Real V1=waterVolume(); + while (abs((totalVolumeofWater-V1)/totalVolumeofWater)>eps){ + slope= (p1-p0)/(V1-V0); + p0=p1; + V0=V1; + p1=p1-slope*(V1-totalVolumeofWater); + if (p1<0) { + cout<< "The requested volume of water is quite big, the simulation will continue at constant suction.:"<< p0 <<endl; + capillaryPressure=p0; + imposePressure=1; + break; + } + solver(p1,switched); + V1=waterVolume(); + capillaryPressure=p1; + } + if (x==1){ + totalVolumeConstant=0; + firstIteration+=1; + } + } + else{ + if ((!totalVolumeConstant && firstIteration==1) && totalVolumeofWater==-1){ + totalVolumeConstant=1; + solver(capillaryPressure,switched); + firstIteration+=1; + totalVolumeConstant=0; + } + else + { + solver(capillaryPressure,switched); + } + + } + } + for (ii= scene->interactions->begin(); ii!=iiEnd ; ++ii) { + CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>((*ii)->phys.get()); + if ((*ii)->isReal() && phys-> computeBridge==true) { + CapillaryPhys1* cundallContactPhysics=NULL; + MindlinCapillaryPhys* mindlinContactPhysics=NULL; + if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys1*>((*ii)->phys.get());//use CapillaryPhys for linear model + else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>((*ii)->phys.get());//use MindlinCapillaryPhys for hertz model + + if ((hertzOn && mindlinContactPhysics->meniscus) || (!hertzOn && cundallContactPhysics->meniscus)) { + if (fusionDetection) {//version with effect of fusion +//BINARY VERSION : if fusionNumber!=0 then no capillary force + short int& fusionNumber = hertzOn?mindlinContactPhysics->fusionNumber:cundallContactPhysics->fusionNumber; + if (binaryFusion) { + if (fusionNumber!=0) { //cerr << "fusion" << endl; + 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.); + } + scene->forces.addForce((*ii)->getId1(), hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap); + scene->forces.addForce((*ii)->getId2(),-(hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap)); + } + } + } + +} +void Law2_ScGeom_CapillaryPhys_Capillarity1::checkFusion() +{ +//Reset fusion numbers + InteractionContainer::iterator ii = scene->interactions->begin(); + InteractionContainer::iterator iiEnd = scene->interactions->end(); + for( ; ii!=iiEnd ; ++ii ) { + if ((*ii)->isReal()) { + if (!hertzOn) static_cast<CapillaryPhys1*>((*ii)->phys.get())->fusionNumber=0; + else static_cast<MindlinCapillaryPhys*>((*ii)->phys.get())->fusionNumber=0; + } + } + + std::list< shared_ptr<Interaction> >::iterator firstMeniscus, lastMeniscus, currentMeniscus; + Real angle1 = -1.0; + Real angle2 = -1.0; + + for ( int i=0; i< bodiesMenisciiList.size(); ++i ) { // i is the index (or id) of the body being tested + CapillaryPhys1* cundallInteractionPhysics1=NULL; + MindlinCapillaryPhys* mindlinInteractionPhysics1=NULL; + CapillaryPhys1* cundallInteractionPhysics2=NULL; + MindlinCapillaryPhys* mindlinInteractionPhysics2=NULL; + if ( !bodiesMenisciiList[i].empty() ) { + lastMeniscus = bodiesMenisciiList[i].end(); + for ( firstMeniscus=bodiesMenisciiList[i].begin(); firstMeniscus!=lastMeniscus; ++firstMeniscus ) { //FOR EACH MENISCUS ON THIS BODY... + currentMeniscus = firstMeniscus; + ++currentMeniscus; + if (!hertzOn) { + cundallInteractionPhysics1 = YADE_CAST<CapillaryPhys1*>((*firstMeniscus)->phys.get()); + if (i == (*firstMeniscus)->getId1()) angle1=cundallInteractionPhysics1->Delta1;//get angle of meniscus1 on body i + else angle1=cundallInteractionPhysics1->Delta2; + } + else { + mindlinInteractionPhysics1 = YADE_CAST<MindlinCapillaryPhys*>((*firstMeniscus)->phys.get()); + if (i == (*firstMeniscus)->getId1()) angle1=mindlinInteractionPhysics1->Delta1;//get angle of meniscus1 on body i + else angle1=mindlinInteractionPhysics1->Delta2; + } + for ( ; currentMeniscus!= lastMeniscus; ++currentMeniscus) { //... CHECK FUSION WITH ALL OTHER MENISCII ON THE BODY + if (!hertzOn) { + cundallInteractionPhysics2 = YADE_CAST<CapillaryPhys1*>((*currentMeniscus)->phys.get()); + if (i == (*currentMeniscus)->getId1()) angle2=cundallInteractionPhysics2->Delta1;//get angle of meniscus2 on body i + else angle2=cundallInteractionPhysics2->Delta2; + } + else { + mindlinInteractionPhysics2 = YADE_CAST<MindlinCapillaryPhys*>((*currentMeniscus)->phys.get()); + if (i == (*currentMeniscus)->getId1()) angle2=mindlinInteractionPhysics2->Delta1;//get angle of meniscus2 on body i + else angle2=mindlinInteractionPhysics2->Delta2; + } + if (angle1==0 || angle2==0) cerr << "THIS SHOULD NOT HAPPEN!!"<< endl; + + + Vector3r normalFirstMeniscus = YADE_CAST<ScGeom*>((*firstMeniscus)->geom.get())->normal; + Vector3r normalCurrentMeniscus = YADE_CAST<ScGeom*>((*currentMeniscus)->geom.get())->normal; + + Real normalDot = 0; + if ((*firstMeniscus)->getId1() == (*currentMeniscus)->getId1() || (*firstMeniscus)->getId2() == (*currentMeniscus)->getId2()) normalDot = normalFirstMeniscus.dot(normalCurrentMeniscus); + else normalDot = - (normalFirstMeniscus.dot(normalCurrentMeniscus)); + + Real normalAngle = 0; + if (normalDot >= 0 ) normalAngle = Mathr::FastInvCos0(normalDot); + else normalAngle = ((Mathr::PI) - Mathr::FastInvCos0(-(normalDot))); + + if ((angle1+angle2)*Mathr::DEG_TO_RAD > normalAngle) { + if (!hertzOn) { + ++(cundallInteractionPhysics1->fusionNumber); //count +1 if 2 meniscii are overlaping + ++(cundallInteractionPhysics2->fusionNumber); + } + else { + ++(mindlinInteractionPhysics1->fusionNumber); + ++(mindlinInteractionPhysics2->fusionNumber); + } + }; + } + } + } + } +} + + +BodiesMenisciiList1::BodiesMenisciiList1(Scene * scene) +{ + initialized=false; + prepare(scene); +} + +bool BodiesMenisciiList1::prepare(Scene * scene) +{ + interactionsOnBody.clear(); + shared_ptr<BodyContainer>& bodies = scene->bodies; + + Body::id_t MaxId = -1; + BodyContainer::iterator bi = bodies->begin(); + BodyContainer::iterator biEnd = bodies->end(); + for( ; bi!=biEnd ; ++bi ) + { + MaxId=max(MaxId, (*bi)->getId()); + } + interactionsOnBody.resize(MaxId+1); + for ( unsigned int i=0; i<interactionsOnBody.size(); ++i ) + { + interactionsOnBody[i].clear(); + } + + InteractionContainer::iterator ii = scene->interactions->begin(); + InteractionContainer::iterator iiEnd = scene->interactions->end(); + for( ; ii!=iiEnd ; ++ii ) { + if ((*ii)->isReal()) { + if (static_cast<CapillaryPhys1*>((*ii)->phys.get())->meniscus) insert(*ii); + } + } + + return initialized=true; +} + +bool BodiesMenisciiList1::insert(const shared_ptr< Interaction >& interaction) +{ + interactionsOnBody[interaction->getId1()].push_back(interaction); + interactionsOnBody[interaction->getId2()].push_back(interaction); + return true; +} + + +bool BodiesMenisciiList1::remove(const shared_ptr< Interaction >& interaction) +{ + interactionsOnBody[interaction->getId1()].remove(interaction); + interactionsOnBody[interaction->getId2()].remove(interaction); + return true; +} + +list< shared_ptr<Interaction> >& BodiesMenisciiList1::operator[] (int index) +{ + return interactionsOnBody[index]; +} + +int BodiesMenisciiList1::size() +{ + return interactionsOnBody.size(); +} + +void BodiesMenisciiList1::display() +{ + list< shared_ptr<Interaction> >::iterator firstMeniscus; + list< shared_ptr<Interaction> >::iterator lastMeniscus; + for ( unsigned int i=0; i<interactionsOnBody.size(); ++i ) + { + if ( !interactionsOnBody[i].empty() ) + { + lastMeniscus = interactionsOnBody[i].end(); + for ( firstMeniscus=interactionsOnBody[i].begin(); firstMeniscus!=lastMeniscus; ++firstMeniscus ) + { + if ( *firstMeniscus ) { + if ( firstMeniscus->get() ) + cerr << "(" << ( *firstMeniscus )->getId1() << ", " << ( *firstMeniscus )->getId2() <<") "; + else cerr << "(void)"; + } + } + cerr << endl; + } + else cerr << "empty" << endl; + } +} + +BodiesMenisciiList1::BodiesMenisciiList1() +{ + initialized=false; +} + +void Law2_ScGeom_CapillaryPhys_Capillarity1::solver(Real suction,bool reset) +{ + + if (!scene) cerr << "scene not defined!"; + shared_ptr<BodyContainer>& bodies = scene->bodies; + if (dtPbased.number_of_vertices ()<1 ) triangulateData(); + if (fusionDetection && !bodiesMenisciiList.initialized) bodiesMenisciiList.prepare(scene); + InteractionContainer::iterator ii = scene->interactions->begin(); + InteractionContainer::iterator iiEnd = scene->interactions->end(); + bool hertzInitialized = false; + Real i=0; + + for (; ii!=iiEnd ; ++ii) { + + i+=1; + CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>((*ii)->phys.get());//////////////////////////////////////////////////////////////////////////////////////////// + + + if ((*ii)->isReal() && phys-> computeBridge==true) { + const shared_ptr<Interaction>& interaction = *ii; + if (!hertzInitialized) {//NOTE: We are assuming that only one type is used in one simulation here + if (CapillaryPhys1::getClassIndexStatic()==interaction->phys->getClassIndex()) hertzOn=false; + else if (MindlinCapillaryPhys::getClassIndexStatic()==interaction->phys->getClassIndex()) hertzOn=true; + else LOG_ERROR("The capillary law is not implemented for interactions using"<<interaction->phys->getClassName()); + } + hertzInitialized = true; + CapillaryPhys1* cundallContactPhysics=NULL; + MindlinCapillaryPhys* mindlinContactPhysics=NULL; + +/// contact physics depends on the contact law, that is used (either linear model or hertz model) + if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys1*>(interaction->phys.get());//use CapillaryPhys for linear model + else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>(interaction->phys.get());//use MindlinCapillaryPhys for hertz model + + unsigned int id1 = interaction->getId1(); + unsigned int id2 = interaction->getId2(); + +/// interaction geometry search (this test is to compute capillarity only between spheres (probably a better way to do that) + int geometryIndex1 = (*bodies)[id1]->shape->getClassIndex(); // !!! + int geometryIndex2 = (*bodies)[id2]->shape->getClassIndex(); + if (!(geometryIndex1 == geometryIndex2)) continue; + +/// definition of interacting objects (not necessarily in contact) + ScGeom* currentContactGeometry = static_cast<ScGeom*>(interaction->geom.get()); + +/// Interacting Grains: +// If you want to define a ratio between YADE sphere size and real sphere size + Real alpha=1; + Real R1 = alpha*std::max(currentContactGeometry->radius2,currentContactGeometry->radius1); + Real R2 =alpha*std::min(currentContactGeometry->radius2,currentContactGeometry->radius1); + + shared_ptr<BodyContainer>& bodies = scene->bodies; + + Real N=bodies->size(); + Real epsilon1,epsilon2; + Real ran1= (N-id1+0.5)/(N+1); + Real ran2= (N-id2+0.5)/(N+1); + epsilon1 = epsilonMean*(2*(ran1-0.5)* disp +1);//addednow + epsilon2 = epsilonMean*(2*(ran2-0.5)* disp +1); +// cout << epsilon1 << "separate" <<epsilon2 <<endl; + R1 = R1-epsilon1*R1; + R2 =R2-epsilon2*R2; + +/// intergranular distance + Real D = alpha*(-(currentContactGeometry->penetrationDepth))+epsilon1*R1+epsilon2*R2; + if ((currentContactGeometry->penetrationDepth>=0)|| D<=0 || createDistantMeniscii) { //||(scene->iter < 1) ) // a simplified way to define meniscii everywhere + if (!hertzOn) { + if (fusionDetection && !cundallContactPhysics->meniscus) bodiesMenisciiList.insert((*ii)); + cundallContactPhysics->meniscus=true; + } else { + if (fusionDetection && !mindlinContactPhysics->meniscus) bodiesMenisciiList.insert((*ii)); + mindlinContactPhysics->meniscus=true; + } + + + } + Real Dinterpol = D/R1; + +/// Suction (Capillary pressure): + if (imposePressure || (!imposePressure && totalVolumeConstant)){ + Real Pinterpol = 0; + if (!hertzOn) Pinterpol = cundallContactPhysics->isBroken ? 0 : suction*R1/liquidTension; + else Pinterpol = mindlinContactPhysics->isBroken ? 0 : suction*R1/liquidTension; + if (!hertzOn) cundallContactPhysics->capillaryPressure = suction; + else mindlinContactPhysics->capillaryPressure = suction; + /// Capillary solution finder: + if ((Pinterpol>=0) && (hertzOn? mindlinContactPhysics->meniscus : cundallContactPhysics->meniscus)) { + MeniscusPhysicalData solution = interpolate1(dtPbased,K::Point_3(R2/R1, Pinterpol, Dinterpol), cundallContactPhysics->m, solutions, reset); +/// capillary adhesion force + Real Finterpol = solution.force; + Vector3r fCap = Finterpol*R1*liquidTension*currentContactGeometry->normal; + if (!hertzOn) cundallContactPhysics->fCap = fCap; + else mindlinContactPhysics->fCap = fCap; +/// meniscus volume +//FIXME: hardcoding numerical constants is bad practice generaly, and it probably reveals a flaw in that case (Bruno) + Real Vinterpol = solution.volume*pow(R1,3); + Real SInterface = solution.surface*pow(R1,2); + if (!hertzOn) { + cundallContactPhysics->vMeniscus = Vinterpol; + cundallContactPhysics->SInterface = SInterface; + if (Vinterpol > 0) cundallContactPhysics->meniscus = true; + else cundallContactPhysics->meniscus = false; + } else { + mindlinContactPhysics->vMeniscus = Vinterpol; + if (Vinterpol > 0) mindlinContactPhysics->meniscus = true; + else mindlinContactPhysics->meniscus = false; + } + if (cundallContactPhysics->meniscus== false) cundallContactPhysics->SInterface=4*3.141592653589793238462643383279502884*(pow(R1,2))+4*3.141592653589793238462643383279502884*(pow(R2,2)); + if (!Vinterpol) { + if ((fusionDetection) || (hertzOn ? mindlinContactPhysics->isBroken : cundallContactPhysics->isBroken)) bodiesMenisciiList.remove((*ii)); + if (D>((interactionDetectionFactor-1)*(currentContactGeometry->radius2+currentContactGeometry->radius1))) scene->interactions->requestErase(interaction); + } +/// wetting angles + if (!hertzOn) { + cundallContactPhysics->Delta1 = max(solution.delta1,solution.delta2); + cundallContactPhysics->Delta2 = min(solution.delta1,solution.delta2); + } else { + mindlinContactPhysics->Delta1 = max(solution.delta1,solution.delta2); + mindlinContactPhysics->Delta2 = min(solution.delta1,solution.delta2); + } + } + + } + else{ + if (cundallContactPhysics->vMeniscus==0 and cundallContactPhysics->capillaryPressure!=0){ + Real Pinterpol = 0; + if (!hertzOn) Pinterpol = cundallContactPhysics->isBroken ? 0 : cundallContactPhysics->capillaryPressure*R1/liquidTension; + else Pinterpol = mindlinContactPhysics->isBroken ? 0 : cundallContactPhysics->capillaryPressure*R1/liquidTension; +// if (!hertzOn) cundallContactPhysics->capillaryPressure = suction; +// else mindlinContactPhysics->capillaryPressure = suction; + /// Capillary solution finder: + if ((Pinterpol>=0) && (hertzOn? mindlinContactPhysics->meniscus : cundallContactPhysics->meniscus)) { + MeniscusPhysicalData solution = interpolate1(dtPbased,K::Point_3(R2/R1, Pinterpol, Dinterpol), cundallContactPhysics->m, solutions, reset); +/// capillary adhesion force + Real Finterpol = solution.force; + Vector3r fCap = Finterpol*R1*liquidTension*currentContactGeometry->normal; + if (!hertzOn) cundallContactPhysics->fCap = fCap; + else mindlinContactPhysics->fCap = fCap; +/// meniscus volume +//FIXME: hardcoding numerical constants is bad practice generaly, and it probably reveals a flaw in that case (Bruno) + Real Vinterpol = solution.volume*pow(R1,3); + Real SInterface = solution.surface*pow(R1,2); + if (!hertzOn) { + cundallContactPhysics->vMeniscus = Vinterpol; + cundallContactPhysics->SInterface = SInterface; + if (Vinterpol > 0) cundallContactPhysics->meniscus = true; + else cundallContactPhysics->meniscus = false; + } else { + mindlinContactPhysics->vMeniscus = Vinterpol; + if (Vinterpol > 0) mindlinContactPhysics->meniscus = true; + else mindlinContactPhysics->meniscus = false; + } + if (cundallContactPhysics->meniscus== false) cundallContactPhysics->SInterface=4*3.141592653589793238462643383279502884*(pow(R1,2))+4*3.141592653589793238462643383279502884*(pow(R2,2)); + if (!Vinterpol) { + if ((fusionDetection) || (hertzOn ? mindlinContactPhysics->isBroken : cundallContactPhysics->isBroken)) bodiesMenisciiList.remove((*ii)); + if (D>((interactionDetectionFactor-1)*(currentContactGeometry->radius2+currentContactGeometry->radius1))) scene->interactions->requestErase(interaction); + } +/// wetting angles + if (!hertzOn) { + cundallContactPhysics->Delta1 = max(solution.delta1,solution.delta2); + cundallContactPhysics->Delta2 = min(solution.delta1,solution.delta2); + } else { + mindlinContactPhysics->Delta1 = max(solution.delta1,solution.delta2); + mindlinContactPhysics->Delta2 = min(solution.delta1,solution.delta2); + } + } + } + else{ + Real Vinterpol = 0; + if (!hertzOn) { Vinterpol = cundallContactPhysics->vMeniscus/pow(R1,3); + } + else Vinterpol = mindlinContactPhysics->vMeniscus; + /// Capillary solution finder: + if ((hertzOn? mindlinContactPhysics->meniscus : cundallContactPhysics->meniscus)) { + MeniscusPhysicalData solution = interpolate2(dtVbased,K::Point_3(R2/R1, Vinterpol, Dinterpol), cundallContactPhysics->m, solutions,reset); +/// capillary adhesion force + Real Finterpol = solution.force; + Vector3r fCap = Finterpol*R1*liquidTension*currentContactGeometry->normal; + if (!hertzOn) cundallContactPhysics->fCap = fCap; + else mindlinContactPhysics->fCap = fCap; +/// suction and interfacial area + + Real Pinterpol = solution.succion*liquidTension/R1; + Real SInterface = solution.surface*pow(R1,2); + if (!hertzOn) { + cundallContactPhysics->capillaryPressure = Pinterpol; + cundallContactPhysics->SInterface = SInterface; + if (Finterpol > 0) cundallContactPhysics->meniscus = true; + else{ + cundallContactPhysics->vMeniscus=0; + cundallContactPhysics->meniscus = false; + } + } else { + mindlinContactPhysics->capillaryPressure = Pinterpol; + if (Finterpol > 0) mindlinContactPhysics->meniscus = true; + else { + cundallContactPhysics->vMeniscus=0; + mindlinContactPhysics->meniscus = false; + } + } + if (!Vinterpol) { + if ((fusionDetection) || (hertzOn ? mindlinContactPhysics->isBroken : cundallContactPhysics->isBroken)) bodiesMenisciiList.remove((*ii)); + if (D>((interactionDetectionFactor-1)*(currentContactGeometry->radius2+currentContactGeometry->radius1))) scene->interactions->requestErase(interaction); + } +/// wetting angles + if (!hertzOn) { + cundallContactPhysics->Delta1 = max(solution.delta1,solution.delta2); + cundallContactPhysics->Delta2 = min(solution.delta1,solution.delta2); + } else { + mindlinContactPhysics->Delta1 = max(solution.delta1,solution.delta2); + mindlinContactPhysics->Delta2 = min(solution.delta1,solution.delta2); + } + } + } + } + + + +///interaction is not real //If the interaction is not real, it should not be in the list + } else if (fusionDetection) bodiesMenisciiList.remove((*ii)); + + + } + if (fusionDetection) checkFusion(); + +} + +#endif //LAW2_SCGEOM_CAPILLARYPHYS_Capillarity1 === added file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.hpp' --- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.hpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.hpp 2015-07-23 17:26:21 +0000 @@ -0,0 +1,123 @@ +// +// C++ Interface: Law2_ScGeom_CapillaryPhys_Capillarity +/************************************************************************* +* Copyright (C) 2006 by luc Scholtes * +* luc.schol...@hmg.inpg.fr * +* * +* This program is free software; it is licensed under the terms of the * +* GNU General Public License v2 or later. See file LICENSE for details. * +*************************************************************************/ +#pragma once + +#include <core/GlobalEngine.hpp> +#include <set> +#include <boost/tuple/tuple.hpp> + +#include <vector> +#include <list> +#include <utility> +#include <pkg/dem/CapillaryPhys1.hpp> +#include<pkg/common/Dispatching.hpp> +#include <string> +#include <iostream> +#include <fstream> + +/** +This law allows one to take into account capillary forces/effects between spheres coming from the presence of interparticular liquid bridges (menisci). +refs: +- (french, lot of documentation) L. Scholtes, PhD thesis -> http://tel.archives-ouvertes.fr/tel-00363961/en/ +- (english, less...) L. Scholtes et al. Micromechanics of granular materials with capillary effects. International Journal of Engineering Science 2009,(47)1, 64-75 + +The law needs ascii files M(r=i) with i=R1/R2 to work (downloaded from https://yade-dem.org/wiki/CapillaryTriaxialTest). They contain a set of results from the resolution of the Laplace-Young equation for different configurations of the interacting geometry and must be placed in the bin directory (where yade exec file is situated) to be taken into account. +The control parameter is the capillary pressure (or suction) Delta_u, defined as the difference between gas and liquid pressure: Delta_u = u_gas - u_liquid +Liquid bridges properties (volume V, extent over interacting grains delta1 and delta2) are computed as a result of Delta_u and the interacting geometry (spheres radii and interparticular distance) + +Rk: - the formulation is valid only for pendular menisci involving two grains (pendular regime). +- an algorithm was developed by B. Chareyre to identify menisci overlaps on each spheres (menisci fusion). +- some assumptions can be made to reduce capillary forces when menisci overlap (binary->F_cap=0 if at least 1 overlap, linear->F_cap=F_cap/numberOfOverlaps) +*/ + +/// !!! This version is deprecated. It should be updated to the new formalism -> ToDo !!! + + + +/// R = ratio(RadiusParticle1 on RadiusParticle2). Here, 10 R values from interpolation files (yade/extra/capillaryFiles), R = 1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 +//const int NB_R_VALUES = 10; + +// class capillarylaw1; // fait appel a la classe def plus bas +class Interaction; + + +///This container class is used to check if meniscii overlap. Wet interactions are put in a series of lists, with one list per body. +class BodiesMenisciiList1 +{ +private: +vector< list< shared_ptr<Interaction> > > interactionsOnBody; + +//shared_ptr<Interaction> empty; + +public: +BodiesMenisciiList1(); +BodiesMenisciiList1(Scene* body); +bool prepare(Scene* scene); +bool insert(const shared_ptr<Interaction>& interaction); +bool remove(const shared_ptr<Interaction>& interaction); +list< shared_ptr<Interaction> >& operator[] (int index); +int size(); +void display(); + + +bool initialized; +}; + +/// This is the constitutive law +class Law2_ScGeom_CapillaryPhys_Capillarity1 : public GlobalEngine +{ +public : + void checkFusion(); + + static DT dtVbased; + static DT dtPbased; + std::vector<MeniscusPhysicalData> solutions; + int switchTriangulation;//to detect switches between P-based and V-based data + + + BodiesMenisciiList1 bodiesMenisciiList; + + void action(); + Real intEnergy(); + Real swInterface(); + Real wnInterface(); + Real waterVolume(); + void solver(Real suction, bool reset); + void triangulateData(); + + YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_CapillaryPhys_Capillarity1,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")) + ((Real,totalVolumeofWater,-1.,,"Value of imposed water volume")) + ((Real,liquidTension,0.073,,"Value of the superficial water tension in N/m")) + ((Real,epsilonMean,0.,," Mean Value of the roughness"))//old value was epsilon + ((Real,disp,0.,," Dispersion from the mean Value of the roughness"))//added now + ((Real,interactionDetectionFactor,1.5,,"defines critical distance for deleting interactions. Must be consistent with the Ig2 value.")) + ((bool,fusionDetection,false,,"If true potential menisci overlaps are checked")) + ((bool,initialized,false,," ")) + ((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")) + ((string,inputFilename,string("capillaryfile.txt"),,"the file with meniscus solutions, used for interpolation.")) + ((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,imposePressure,true,," If True, suction is imposed and is constant if not Volume is imposed-Undrained test")) + ((bool,totalVolumeConstant,true,," in undrained test there are 2 options, If True, the total volume of water is imposed,if false the volume of each meniscus is kept constant: in this case capillary pressure can be imposed for initial distribution of meniscus or it is the total volume that can be imposed initially")) + ,switchTriangulation=-1, + .def("intEnergy",&Law2_ScGeom_CapillaryPhys_Capillarity1::intEnergy,"define the energy of interfaces in unsaturated pendular state") +// .def("sninterface",&Law2_ScGeom_CapillaryPhys_Capillarity1::sninterface,"define the amount of solid-non-wetting interfaces in unsaturated pendular state") + .def("swInterface",&Law2_ScGeom_CapillaryPhys_Capillarity1::swInterface,"define the amount of solid-wetting interfaces in unsaturated pendular state") + .def("wnInterface",&Law2_ScGeom_CapillaryPhys_Capillarity1::wnInterface,"define the amount of wetting-non-wetiing interfaces in unsaturated pendular state") + .def("waterVolume",&Law2_ScGeom_CapillaryPhys_Capillarity1::waterVolume,"return the total value of water in the sample") + ); +}; + + + +REGISTER_SERIALIZABLE(Law2_ScGeom_CapillaryPhys_Capillarity1); +
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp