yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10631
[Branch ~yade-pkg/yade/git-trunk] Rev 3866: Move capillary functions out of ViscEl-classes.
------------------------------------------------------------
revno: 3866
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Tue 2014-04-01 16:45:46 +0200
message:
Move capillary functions out of ViscEl-classes.
Discussed here:
http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg09786.html
Thanks Bruno and Klaus.
modified:
pkg/dem/ViscoelasticCapillarPM.cpp
pkg/dem/ViscoelasticCapillarPM.hpp
pkg/dem/ViscoelasticPM.cpp
pkg/dem/ViscoelasticPM.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 2014-04-01 14:45:46 +0000
+++ pkg/dem/ViscoelasticCapillarPM.cpp 2014-04-01 14:45:46 +0000
@@ -5,10 +5,139 @@
#include<yade/core/Scene.hpp>
#include<yade/pkg/common/Sphere.hpp>
-Real Law2_ScGeom_ViscElPhys_Basic::calculateCapillarForce(const ScGeom& geom, ViscElPhys& phys) {
+YADE_PLUGIN((ViscElCapMat)(ViscElCapPhys)(Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys)(Law2_ScGeom_ViscElCapPhys_Basic));
+
+/* ViscElCapMat */
+ViscElCapMat::~ViscElCapMat(){}
+
+/* ViscElCapPhys */
+ViscElCapPhys::~ViscElCapPhys(){}
+
+/* Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys */
+void Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) {
+ // no updates of an existing contact
+ if(interaction->phys) return;
+
+ ViscElPhys* physT = Calculate_ViscElMat_ViscElMat_ViscElPhys(b1, b2, interaction);
+
+ shared_ptr<ViscElCapPhys> phys (new ViscElCapPhys());
+ phys->kn = physT->kn;
+ phys->ks = physT->ks;
+ phys->cn = physT->cn;
+ phys->cs = physT->cs;
+ phys->tangensOfFrictionAngle = physT->tangensOfFrictionAngle;
+ phys->shearForce = physT->shearForce;
+ phys->mRtype = physT->mRtype;
+
+ ViscElCapMat* mat1 = static_cast<ViscElCapMat*>(b1.get());
+ ViscElCapMat* mat2 = static_cast<ViscElCapMat*>(b2.get());
+
+ if (mat1->Capillar and mat2->Capillar) {
+ if (mat1->Vb == mat2->Vb) {
+ phys->Vb = mat1->Vb;
+ } else {
+ throw runtime_error("Vb should be equal for both particles!.");
+ }
+
+ if (mat1->gamma == mat2->gamma) {
+ phys->gamma = mat1->gamma;
+ } else {
+ throw runtime_error("Gamma should be equal for both particles!.");
+ }
+
+ if (mat1->theta == mat2->theta) {
+ phys->theta = (mat1->theta*M_PI/180.0);
+ } 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;
+ else if (mat1->CapillarType == "Willett_analytic") phys->CapillarType = Willett_analytic;
+ else if (mat1->CapillarType == "Weigert") phys->CapillarType = Weigert;
+ else if (mat1->CapillarType == "Rabinovich") phys->CapillarType = Rabinovich;
+ else if (mat1->CapillarType == "Lambert") phys->CapillarType = Lambert;
+ else phys->CapillarType = None_Capillar;
+ } else {
+ throw runtime_error("CapillarType should be equal for both particles!.");
+ }
+ phys->Capillar=true;
+ }
+
+ interaction->phys = phys;
+}
+
+/* Law2_ScGeom_ViscElCapPhys_Basic */
+void Law2_ScGeom_ViscElCapPhys_Basic::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I) {
+ Vector3r force = Vector3r::Zero();
+
+ const int id1 = I->getId1();
+ const int id2 = I->getId2();
+
+ const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
+ Scene* scene=Omega::instance().getScene().get();
+ ViscElCapPhys& phys=*static_cast<ViscElCapPhys*>(_phys.get());
+ const BodyContainer& bodies = *scene->bodies;
+
+ /*
+ * This part for implementation of the capillar model.
+ * All main equations are in calculateCapillarForce function.
+ * There is only the determination of critical distance between spheres,
+ * after that the liquid bridge will be broken.
+ */
+
+ if (not(phys.liqBridgeCreated) and phys.Capillar) {
+ phys.liqBridgeCreated = true;
+ Sphere* s1=dynamic_cast<Sphere*>(bodies[id1]->shape.get());
+ Sphere* s2=dynamic_cast<Sphere*>(bodies[id2]->shape.get());
+ if (s1 and s2) {
+ phys.R = 2 * s1->radius * s2->radius / (s1->radius + s2->radius);
+ } else if (s1 and not(s2)) {
+ phys.R = s1->radius;
+ } else {
+ phys.R = s2->radius;
+ }
+
+ const Real Vstar = phys.Vb/(phys.R*phys.R*phys.R);
+ const Real Sstar = (1+0.5*phys.theta)*(pow(Vstar,1/3.0) + 0.1*pow(Vstar,2.0/3.0)); // [Willett2000], equation (15), use the full-length e.g 2*Sc
+
+ phys.sCrit = Sstar*phys.R;
+ }
+
+ if (geom.penetrationDepth<0) {
+ if (phys.liqBridgeCreated and -geom.penetrationDepth<phys.sCrit and phys.Capillar) {
+ phys.normalForce = -calculateCapillarForce(geom, phys)*geom.normal;
+ if (I->isActive) {
+ addForce (id1,-phys.normalForce,scene);
+ addForce (id2, phys.normalForce,scene);
+ };
+ return;
+ } else {
+ scene->interactions->requestErase(I);
+ return;
+ };
+ };
+
+ if (I->isActive) {
+
+ Vector3r torque1 = Vector3r::Zero();
+ Vector3r torque2 = Vector3r::Zero();
+
+ computeForceTorqueViscEl(_geom, _phys, I, force, torque1, torque2);
+
+ addForce (id1,-force,scene);
+ addForce (id2, force,scene);
+ addTorque(id1, torque1,scene);
+ addTorque(id2, torque2,scene);
+ }
+}
+
+Real Law2_ScGeom_ViscElCapPhys_Basic::calculateCapillarForce(const ScGeom& geom, ViscElCapPhys& phys) {
Real fC = 0.0;
- /* Capillar
+ /*
+ * Capillar
* Some equations have constants, which can be calculated only once per contact.
* No need to recalculate them at each step.
* It needs to be fixed.
=== modified file 'pkg/dem/ViscoelasticCapillarPM.hpp'
--- pkg/dem/ViscoelasticCapillarPM.hpp 2014-04-01 14:45:46 +0000
+++ pkg/dem/ViscoelasticCapillarPM.hpp 2014-04-01 14:45:46 +0000
@@ -1,1 +1,59 @@
#include"ViscoelasticPM.hpp"
+
+class ViscElCapMat : public ViscElMat {
+ public:
+ virtual ~ViscElCapMat();
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(ViscElCapMat,ViscElMat,"Material for extended viscoelastic model of contact with capillary parameters.",
+ ((bool,Capillar,false,,"True, if capillar forces need to be added."))
+ ((Real,Vb,NaN,,"Liquid bridge volume [m^3]"))
+ ((Real,gamma,NaN,,"Surface tension [N/m]"))
+ ((Real,theta,NaN,,"Contact angle [°]"))
+ ((std::string,CapillarType,"",,"Different types of capillar interaction: Willett_numeric, Willett_analytic [Willett2000]_ , Weigert [Weigert1999]_ , Rabinovich [Rabinov2005]_ , Lambert (simplified, corrected Rabinovich model) [Lambert2008]_ ")),
+ createIndex();
+ );
+ REGISTER_CLASS_INDEX(ViscElCapMat,ViscElMat);
+};
+REGISTER_SERIALIZABLE(ViscElCapMat);
+
+/// Interaction physics
+enum CapType {None_Capillar, Willett_numeric, Willett_analytic, Weigert, Rabinovich, Lambert};
+class ViscElCapPhys : public ViscElPhys{
+ public:
+ virtual ~ViscElCapPhys();
+ Real R;
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(ViscElCapPhys,ViscElPhys,"IPhys created from :yref:`ViscElCapMat`, for use with :yref:`Law2_ScGeom_ViscElCapPhys_Basic`.",
+ ((bool,Capillar,false,,"True, if capillar forces need to be added."))
+ ((bool,liqBridgeCreated,false,,"Whether liquid bridge was created, only after a normal contact of spheres"))
+ ((Real,sCrit,false,,"Critical bridge length [m]"))
+ ((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")),
+ createIndex();
+ )
+ REGISTER_CLASS_INDEX(ViscElCapPhys,ViscElPhys);
+};
+REGISTER_SERIALIZABLE(ViscElCapPhys);
+
+/// Convert material to interaction physics.
+class Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys: public IPhysFunctor {
+ public :
+ virtual void go(const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+ YADE_CLASS_BASE_DOC(Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys,IPhysFunctor,"Convert 2 instances of :yref:`ViscElCapMat` to :yref:`ViscElCapPhys` using the rule of consecutive connection.");
+ FUNCTOR2D(ViscElCapMat,ViscElCapMat);
+};
+REGISTER_SERIALIZABLE(Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys);
+
+/// Constitutive law
+class Law2_ScGeom_ViscElCapPhys_Basic: public LawFunctor {
+ public :
+ virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ public :
+ Real calculateCapillarForce(const ScGeom& geom, ViscElCapPhys& phys);
+ FUNCTOR2D(ScGeom,ViscElCapPhys);
+ YADE_CLASS_BASE_DOC(Law2_ScGeom_ViscElCapPhys_Basic,LawFunctor,"Extended version of Linear viscoelastic model with capillary parameters.");
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom_ViscElCapPhys_Basic);
=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp 2014-04-01 14:45:46 +0000
+++ pkg/dem/ViscoelasticPM.cpp 2014-04-01 14:45:46 +0000
@@ -14,111 +14,11 @@
/* ViscElPhys */
ViscElPhys::~ViscElPhys(){}
-/* Contact parameter calculation function */
-Real Ip2_ViscElMat_ViscElMat_ViscElPhys::contactParameterCalculation(const Real& l1, const Real& l2, const bool& massMultiply){
- if (massMultiply) {
- // If one of paramaters > 0. we DO NOT return 0
- Real a = (l1?1/l1:0) + (l2?1/l2:0);
- if (a) return 1/a;
- else return 0;
- } else {
- // If one of paramaters > 0, we return 0
- return (l1>0 or l2>0)?l1*l2/(l1+l2):0;
- }
-}
-
/* Ip2_ViscElMat_ViscElMat_ViscElPhys */
void Ip2_ViscElMat_ViscElMat_ViscElPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) {
-
// no updates of an existing contact
if(interaction->phys) return;
- ViscElMat* mat1 = static_cast<ViscElMat*>(b1.get());
- ViscElMat* mat2 = static_cast<ViscElMat*>(b2.get());
- Real mass1 = 1.0;
- Real mass2 = 1.0;
-
- if (mat1->massMultiply and mat2->massMultiply) {
- mass1 = Body::byId(interaction->getId1())->state->mass;
- mass2 = Body::byId(interaction->getId2())->state->mass;
- if (mass1 == 0.0 and mass2 > 0.0) {
- mass1 = mass2;
- } else if (mass2 == 0.0 and mass1 > 0.0) {
- mass2 = mass1;
- }
- }
-
- GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
- Real R1=sphCont->refR1>0?sphCont->refR1:sphCont->refR2;
- Real R2=sphCont->refR2>0?sphCont->refR2:sphCont->refR1;
-
- const Real kn1 = isnan(mat1->kn)?2*mat1->young*R1:mat1->kn*mass1;
- const Real cn1 = mat1->cn*mass1;
- const Real ks1 = isnan(mat1->ks)?kn1*mat1->poisson:mat1->ks*mass1;
- const Real cs1 = mat1->cs*mass1;
-
- const Real mR1 = mat1->mR; const Real mR2 = mat2->mR;
- const int mRtype1 = mat1->mRtype; const int mRtype2 = mat2->mRtype;
-
- const Real kn2 = isnan(mat2->kn)?2*mat2->young*R2:mat2->kn*mass2;
- const Real cn2 = mat2->cn*mass2;
- const Real ks2 = isnan(mat2->ks)?kn2*mat2->poisson:mat2->ks*mass2;
- const Real cs2 = mat2->cs*mass2;
-
- ViscElPhys* phys = new ViscElPhys();
-
- phys->kn = contactParameterCalculation(kn1,kn2, mat1->massMultiply&&mat2->massMultiply);
- phys->ks = contactParameterCalculation(ks1,ks2, mat1->massMultiply&&mat2->massMultiply);
- phys->cn = contactParameterCalculation(cn1,cn2, mat1->massMultiply&&mat2->massMultiply);
- phys->cs = contactParameterCalculation(cs1,cs2, mat1->massMultiply&&mat2->massMultiply);
-
- if ((mR1>0) or (mR2>0)) {
- phys->mR = 2.0/( ((mR1>0)?1/mR1:0) + ((mR2>0)?1/mR2:0) );
- } else {
- phys->mR = 0;
- }
-
- phys->tangensOfFrictionAngle = std::tan(std::min(mat1->frictionAngle, mat2->frictionAngle));
- phys->shearForce = Vector3r(0,0,0);
-
- if ((mRtype1 != mRtype2) or (mRtype1>2) or (mRtype2>2) or (mRtype1<1) or (mRtype2<1) ) {
- throw runtime_error("mRtype should be equal for both materials and have the values 1 or 2!");
- } else {
- phys->mRtype = mRtype1;
- }
-
- if (mat1->Capillar and mat2->Capillar) {
- if (mat1->Vb == mat2->Vb) {
- phys->Vb = mat1->Vb;
- } else {
- throw runtime_error("Vb should be equal for both particles!.");
- }
-
- if (mat1->gamma == mat2->gamma) {
- phys->gamma = mat1->gamma;
- } else {
- throw runtime_error("Gamma should be equal for both particles!.");
- }
-
- if (mat1->theta == mat2->theta) {
- phys->theta = (mat1->theta*M_PI/180.0);
- } 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;
- else if (mat1->CapillarType == "Willett_analytic") phys->CapillarType = Willett_analytic;
- else if (mat1->CapillarType == "Weigert") phys->CapillarType = Weigert;
- else if (mat1->CapillarType == "Rabinovich") phys->CapillarType = Rabinovich;
- else if (mat1->CapillarType == "Lambert") phys->CapillarType = Lambert;
- else phys->CapillarType = None_Capillar;
- } else {
- throw runtime_error("CapillarType should be equal for both particles!.");
- }
- phys->Capillar=true;
- }
-
+ ViscElPhys * phys = Calculate_ViscElMat_ViscElMat_ViscElPhys(b1, b2, interaction);
interaction->phys = shared_ptr<ViscElPhys>(phys);
}
@@ -127,7 +27,7 @@
Vector3r force = Vector3r::Zero();
Vector3r torque1 = Vector3r::Zero();
Vector3r torque2 = Vector3r::Zero();
- computeForceTorque(_geom, _phys, I, force, torque1, torque2);
+ computeForceTorqueViscEl(_geom, _phys, I, force, torque1, torque2);
if (I->isActive) {
const int id1 = I->getId1();
const int id2 = I->getId2();
@@ -139,57 +39,19 @@
}
}
-void Law2_ScGeom_ViscElPhys_Basic::computeForceTorque(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2) {
+void computeForceTorqueViscEl(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2) {
const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
+ Scene* scene=Omega::instance().getScene().get();
ViscElPhys& phys=*static_cast<ViscElPhys*>(_phys.get());
const int id1 = I->getId1();
const int id2 = I->getId2();
-
- if (geom.penetrationDepth<0) {
- if (phys.liqBridgeCreated and -geom.penetrationDepth<phys.sCrit and phys.Capillar) {
- phys.normalForce = -calculateCapillarForce(geom, phys)*geom.normal;
- if (I->isActive) {
- addForce (id1,-phys.normalForce,scene);
- addForce (id2, phys.normalForce,scene);
- };
- return;
- } else {
- scene->interactions->requestErase(I);
- return;
- };
- };
const BodyContainer& bodies = *scene->bodies;
const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
const State& de2 = *static_cast<State*>(bodies[id2]->state.get());
- /*
- * This part for implementation of the capillar model.
- * All main equations are in calculateCapillarForce function.
- * There is only the determination of critical distance between spheres,
- * after that the liquid bridge will be broken.
- */
-
- if (not(phys.liqBridgeCreated) and phys.Capillar) {
- phys.liqBridgeCreated = true;
- Sphere* s1=dynamic_cast<Sphere*>(bodies[id1]->shape.get());
- Sphere* s2=dynamic_cast<Sphere*>(bodies[id2]->shape.get());
- if (s1 and s2) {
- phys.R = 2 * s1->radius * s2->radius / (s1->radius + s2->radius);
- } else if (s1 and not(s2)) {
- phys.R = s1->radius;
- } else {
- phys.R = s2->radius;
- }
-
- const Real Vstar = phys.Vb/(phys.R*phys.R*phys.R);
- const Real Sstar = (1+0.5*phys.theta)*(pow(Vstar,1/3.0) + 0.1*pow(Vstar,2.0/3.0)); // [Willett2000], equation (15), use the full-length e.g 2*Sc
-
- phys.sCrit = Sstar*phys.R;
- }
-
Vector3r& shearForce = phys.shearForce;
if (I->isFresh(scene)) shearForce=Vector3r(0,0,0);
const Real& dt = scene->dt;
@@ -253,3 +115,76 @@
torque1 = -c1x.cross(force)+momentResistance;
torque2 = c2x.cross(force)-momentResistance;
}
+
+ViscElPhys* Calculate_ViscElMat_ViscElMat_ViscElPhys(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) {
+ ViscElPhys* phys = new ViscElPhys();
+
+ ViscElMat* mat1 = static_cast<ViscElMat*>(b1.get());
+ ViscElMat* mat2 = static_cast<ViscElMat*>(b2.get());
+ Real mass1 = 1.0;
+ Real mass2 = 1.0;
+
+ if (mat1->massMultiply and mat2->massMultiply) {
+ mass1 = Body::byId(interaction->getId1())->state->mass;
+ mass2 = Body::byId(interaction->getId2())->state->mass;
+ if (mass1 == 0.0 and mass2 > 0.0) {
+ mass1 = mass2;
+ } else if (mass2 == 0.0 and mass1 > 0.0) {
+ mass2 = mass1;
+ }
+ }
+
+ GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
+ Real R1=sphCont->refR1>0?sphCont->refR1:sphCont->refR2;
+ Real R2=sphCont->refR2>0?sphCont->refR2:sphCont->refR1;
+
+ const Real kn1 = isnan(mat1->kn)?2*mat1->young*R1:mat1->kn*mass1;
+ const Real cn1 = mat1->cn*mass1;
+ const Real ks1 = isnan(mat1->ks)?kn1*mat1->poisson:mat1->ks*mass1;
+ const Real cs1 = mat1->cs*mass1;
+
+ const Real mR1 = mat1->mR; const Real mR2 = mat2->mR;
+ const int mRtype1 = mat1->mRtype; const int mRtype2 = mat2->mRtype;
+
+ const Real kn2 = isnan(mat2->kn)?2*mat2->young*R2:mat2->kn*mass2;
+ const Real cn2 = mat2->cn*mass2;
+ const Real ks2 = isnan(mat2->ks)?kn2*mat2->poisson:mat2->ks*mass2;
+ const Real cs2 = mat2->cs*mass2;
+
+
+ phys->kn = contactParameterCalculation(kn1,kn2, mat1->massMultiply&&mat2->massMultiply);
+ phys->ks = contactParameterCalculation(ks1,ks2, mat1->massMultiply&&mat2->massMultiply);
+ phys->cn = contactParameterCalculation(cn1,cn2, mat1->massMultiply&&mat2->massMultiply);
+ phys->cs = contactParameterCalculation(cs1,cs2, mat1->massMultiply&&mat2->massMultiply);
+
+ if ((mR1>0) or (mR2>0)) {
+ phys->mR = 2.0/( ((mR1>0)?1/mR1:0) + ((mR2>0)?1/mR2:0) );
+ } else {
+ phys->mR = 0;
+ }
+
+ phys->tangensOfFrictionAngle = std::tan(std::min(mat1->frictionAngle, mat2->frictionAngle));
+ phys->shearForce = Vector3r(0,0,0);
+
+ if ((mRtype1 != mRtype2) or (mRtype1>2) or (mRtype2>2) or (mRtype1<1) or (mRtype2<1) ) {
+ throw runtime_error("mRtype should be equal for both materials and have the values 1 or 2!");
+ } else {
+ phys->mRtype = mRtype1;
+ }
+
+ return phys;
+}
+
+/* Contact parameter calculation function */
+Real contactParameterCalculation(const Real& l1, const Real& l2, const bool& massMultiply){
+ if (massMultiply) {
+ // If one of paramaters > 0. we DO NOT return 0
+ Real a = (l1?1/l1:0) + (l2?1/l2:0);
+ if (a) return 1/a;
+ else return 0;
+ } else {
+ // If one of paramaters > 0, we return 0
+ return (l1>0 or l2>0)?l1*l2/(l1+l2):0;
+ }
+}
+
=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp 2014-04-01 14:45:46 +0000
+++ pkg/dem/ViscoelasticPM.hpp 2014-04-01 14:45:46 +0000
@@ -37,7 +37,6 @@
REGISTER_SERIALIZABLE(ViscElMat);
/// Interaction physics
-enum CapType {None_Capillar, Willett_numeric, Willett_analytic, Weigert, Rabinovich, Lambert};
class ViscElPhys : public FrictPhys{
public:
virtual ~ViscElPhys();
@@ -46,14 +45,7 @@
((Real,cn,NaN,,"Normal viscous constant"))
((Real,cs,NaN,,"Shear viscous constant"))
((Real,mR,0.0,,"Rolling resistance, see [Zhou1999536]_."))
- ((unsigned int,mRtype,1,,"Rolling resistance type, see [Zhou1999536]_. mRtype=1 - equation (3) in [Zhou1999536]_; mRtype=2 - equation (4) in [Zhou1999536]_"))
- ((bool,Capillar,false,,"True, if capillar forces need to be added."))
- ((bool,liqBridgeCreated,false,,"Whether liquid bridge was created, only after a normal contact of spheres"))
- ((Real,sCrit,false,,"Critical bridge length [m]"))
- ((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")),
+ ((unsigned int,mRtype,1,,"Rolling resistance type, see [Zhou1999536]_. mRtype=1 - equation (3) in [Zhou1999536]_; mRtype=2 - equation (4) in [Zhou1999536]_")),
createIndex();
)
REGISTER_CLASS_INDEX(ViscElPhys,FrictPhys);
@@ -67,8 +59,6 @@
virtual void go(const shared_ptr<Material>& b1,
const shared_ptr<Material>& b2,
const shared_ptr<Interaction>& interaction);
- private :
- Real contactParameterCalculation(const Real& l1,const Real& l2, const bool& massMultiply);
YADE_CLASS_BASE_DOC(Ip2_ViscElMat_ViscElMat_ViscElPhys,IPhysFunctor,"Convert 2 instances of :yref:`ViscElMat` to :yref:`ViscElPhys` using the rule of consecutive connection.");
FUNCTOR2D(ViscElMat,ViscElMat);
@@ -80,9 +70,6 @@
class Law2_ScGeom_ViscElPhys_Basic: public LawFunctor {
public :
virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
- private:
- Real calculateCapillarForce(const ScGeom& geom, ViscElPhys& phys);
- void computeForceTorque(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2);
public :
FUNCTOR2D(ScGeom,ViscElPhys);
YADE_CLASS_BASE_DOC(Law2_ScGeom_ViscElPhys_Basic,LawFunctor,"Linear viscoelastic model operating on :yref:`ScGeom` and :yref:`ViscElPhys`. The model is mostly based on the paper for For details see Pournin [Pournin2001]_ .");
@@ -90,3 +77,6 @@
};
REGISTER_SERIALIZABLE(Law2_ScGeom_ViscElPhys_Basic);
+Real contactParameterCalculation(const Real& l1,const Real& l2, const bool& massMultiply);
+ViscElPhys* Calculate_ViscElMat_ViscElMat_ViscElPhys(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction);
+void computeForceTorqueViscEl(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2);
