yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #12843
[Branch ~yade-pkg/yade/git-trunk] Rev 3934: Including interface orientation data in CapillaryPhys code
------------------------------------------------------------
revno: 3934
committer: jduriez <jerome.duriez@xxxxxxxxxxx>
timestamp: Wed 2016-09-21 10:20:15 -0600
message:
Including interface orientation data in CapillaryPhys code
modified:
pkg/dem/CapillaryPhys.hpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.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/CapillaryPhys.hpp'
--- pkg/dem/CapillaryPhys.hpp 2015-05-22 05:46:49 +0000
+++ pkg/dem/CapillaryPhys.hpp 2016-09-21 16:20:15 +0000
@@ -26,6 +26,8 @@
((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 produced by the presence of the meniscus. This is the force acting on particle #2"))
((short int,fusionNumber,0.,,"Indicates the number of meniscii that overlap with this one"))
+ ((Real,nn11,0.,,":math:`\\iint_A n_1 n_1 \\, dS = \\iint_A n_2 n_2 \\, dS`, $A$ being the liquid-gas surface of the meniscus, $\\vec n$ the associated normal, and $(1,2,3)$ a local basis with $3$ the meniscus orientation (:yref:`ScGeom.normal`). NB: $A$ = 2 :yref:`nn11<CapillaryPhys.nn11>` + :yref:`nn33<CapillaryPhys.nn33>`."))
+ ((Real,nn33,0.,,":math:`\\iint_A n_3 n_3 \\, dS`, $A$ being the liquid-gas surface of the meniscus, $\\vec n$ the associated normal, and $(1,2,3)$ a local basis with $3$ the meniscus orientation (:yref:`ScGeom.normal`). NB: $A$ = 2 :yref:`nn11<CapillaryPhys.nn11>` + :yref:`nn33<CapillaryPhys.nn33>`."))
,,
createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
,
=== modified file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp 2016-07-25 19:26:13 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp 2016-09-21 16:20:15 +0000
@@ -21,16 +21,16 @@
void Law2_ScGeom_CapillaryPhys_Capillarity::postLoad(Law2_ScGeom_CapillaryPhys_Capillarity&){
capillary = shared_ptr<capillarylaw>(new capillarylaw);
- capillary->fill("M(r=1)");
- capillary->fill("M(r=1.1)");
- capillary->fill("M(r=1.25)");
- capillary->fill("M(r=1.5)");
- capillary->fill("M(r=1.75)");
- capillary->fill("M(r=2)");
- capillary->fill("M(r=3)");
- capillary->fill("M(r=4)");
- capillary->fill("M(r=5)");
- capillary->fill("M(r=10)");
+ capillary->fill(("M(r=1)"+suffCapFiles).c_str());
+ capillary->fill(("M(r=1.1)"+suffCapFiles).c_str());
+ capillary->fill(("M(r=1.25)"+suffCapFiles).c_str());
+ capillary->fill(("M(r=1.5)"+suffCapFiles).c_str());
+ capillary->fill(("M(r=1.75)"+suffCapFiles).c_str());
+ capillary->fill(("M(r=2)"+suffCapFiles).c_str());
+ capillary->fill(("M(r=3)"+suffCapFiles).c_str());
+ capillary->fill(("M(r=4)"+suffCapFiles).c_str());
+ capillary->fill(("M(r=5)"+suffCapFiles).c_str());
+ capillary->fill(("M(r=10)"+suffCapFiles).c_str());
}
@@ -40,6 +40,8 @@
F = 0;
delta1 = 0;
delta2 = 0;
+ nn11 = 0;
+ nn33 = 0;
}
MeniscusParameters::MeniscusParameters(const MeniscusParameters &source)
@@ -48,6 +50,8 @@
F = source.F;
delta1 = source.delta1;
delta2 = source.delta2;
+ nn11 = source.nn11;
+ nn33 = source.nn33;
}
MeniscusParameters::~MeniscusParameters()
@@ -162,12 +166,17 @@
/// wetting angles
if (!hertzOn) {
- cundallContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
+ cundallContactPhysics->Delta1 = min(max(solution.delta1,solution.delta2),90.0); // undesired values greater than 90 degrees (by few degrees) may be present in the capillary files for high r (~ 10) and very low suction and distance
cundallContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
} else {
- mindlinContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
+ mindlinContactPhysics->Delta1 = min(max(solution.delta1,solution.delta2),90.0);
mindlinContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
}
+ // nn11 and nn33
+ if (!hertzOn) {
+ cundallContactPhysics->nn11 = pow(R2/alpha,2) * solution.nn11;
+ cundallContactPhysics->nn33 = pow(R2/alpha,2) * solution.nn33;
+ }
}
///interaction is not real //If the interaction is not real, it should not be in the list
} else if (fusionDetection) bodiesMenisciiList.remove(interaction);
@@ -286,7 +295,7 @@
}
MeniscusParameters capillarylaw::interpolate(Real R1, Real R2, Real D, Real P, int* index)
-{ //cerr << "interpolate" << endl;
+{
if (R1 > R2) {
Real R3 = R1;
R1 = R2;
@@ -320,6 +329,8 @@
result.F = result_inf.F*(1-r) + r*result_sup.F;
result.delta1 = result_inf.delta1*(1-r) + r*result_sup.delta1;
result.delta2 = result_inf.delta2*(1-r) + r*result_sup.delta2;
+ result.nn11 = result_inf.nn11*(1-r) + r*result_sup.nn11;
+ result.nn33 = result_inf.nn33*(1-r) + r*result_sup.nn33;
i=NB_R_VALUES;
//cerr << "i = " << i << endl;
@@ -329,7 +340,6 @@
{
result = data_complete[i].Interpolate2(D,P, index[0], index[1]);
i=NB_R_VALUES;
- //cerr << "i = " << i << endl;
}
}
return result;
@@ -343,7 +353,6 @@
{
ifstream file (filename);
file >> R;
- //cerr << "r = " << R << endl;
int n_D; //number of D values
file >> n_D;
@@ -373,7 +382,7 @@
MeniscusParameters result_inf;
MeniscusParameters result_sup;
- for ( unsigned int i=0; i < full_data.size(); ++i)
+ for ( unsigned int i=0; i < full_data.size(); ++i) // loop over D values
{
if (full_data[i].D > D ) // ok si D rang�s ds l'ordre croissant
@@ -387,6 +396,8 @@
result.F = result_inf.F*(1-rD) + rD*result_sup.F;
result.delta1 = result_inf.delta1*(1-rD) + rD*result_sup.delta1;
result.delta2 = result_inf.delta2*(1-rD) + rD*result_sup.delta2;
+ result.nn11 = result_inf.nn11*(1-rD) + rD*result_sup.nn11;
+ result.nn33 = result_inf.nn33*(1-rD) + rD*result_sup.nn33;
i = full_data.size();
}
@@ -410,19 +421,20 @@
Real x;
int n_lines; //pb: n_lines is real!!!
file >> n_lines;
- //cout << n_lines << endl;
file.ignore(200, '\n'); // saute les caract�res (200 au maximum) jusque au caract�re \n (fin de ligne)*_
if (n_lines!=0)
for (; i<n_lines; ++i) {
data.push_back(vector<Real> ());
- for (int j=0; j < 6; ++j) // [D,P,V,F,delta1,delta2]
+ for (int j=0; j < 8; ++j) // [D,P,V,F,delta1,delta2,nn11,nn33]
{
file >> x;
data[i].push_back(x);
}
}
+ else
+ LOG_ERROR("Problem regarding the capillary file structure (e.g. n_D is not consistent with the actual data), and/or mismatch with the expected structure by the code ! Will segfault");
D = data[i-1][0];
}
@@ -441,17 +453,23 @@
Real Vinf=data[index-1][2];
Real Delta1inf=data[index-1][4];
Real Delta2inf=data[index-1][5];
+ Real nn11inf = data[index-1][6];
+ Real nn33inf = data[index-1][7];
Real Psup=data[index][1];
Real Fsup=data[index][3];
Real Vsup=data[index][2];
Real Delta1sup=data[index][4];
Real Delta2sup=data[index][5];
+ Real nn11sup = data[index][6];
+ Real nn33sup = data[index][7];
result.V = Vinf+((Vsup-Vinf)/(Psup-Pinf))*(P-Pinf);
result.F = Finf+((Fsup-Finf)/(Psup-Pinf))*(P-Pinf);
result.delta1 = Delta1inf+((Delta1sup-Delta1inf)/(Psup-Pinf))*(P-Pinf);
result.delta2 = Delta2inf+((Delta2sup-Delta2inf)/(Psup-Pinf))*(P-Pinf);
+ result.nn11 = nn11inf + (nn11sup - nn11inf) / (Psup-Pinf) * (P-Pinf);
+ result.nn33 = nn33inf + (nn33sup - nn33inf) / (Psup-Pinf) * (P-Pinf);
return result;
}
@@ -459,37 +477,45 @@
//compteur2+=1;
for (int k=1; k < dataSize; ++k) // Length(data) ??
- { //cerr << "k = " << k << endl;
+ {
if ( data[k][1] > P) // OK si P rangés ds l'ordre croissant
- { //cerr << "if" << endl;
+ {
Real Pinf=data[k-1][1];
Real Finf=data[k-1][3];
Real Vinf=data[k-1][2];
Real Delta1inf=data[k-1][4];
Real Delta2inf=data[k-1][5];
+ Real nn11inf = data[k-1][6];
+ Real nn33inf = data[k-1][7];
Real Psup=data[k][1];
Real Fsup=data[k][3];
Real Vsup=data[k][2];
Real Delta1sup=data[k][4];
Real Delta2sup=data[k][5];
+ Real nn11sup = data[k][6];
+ Real nn33sup = data[k][7];
result.V = Vinf+((Vsup-Vinf)/(Psup-Pinf))*(P-Pinf);
result.F = Finf+((Fsup-Finf)/(Psup-Pinf))*(P-Pinf);
result.delta1 = Delta1inf+((Delta1sup-Delta1inf)/(Psup-Pinf))*(P-Pinf);
result.delta2 = Delta2inf+((Delta2sup-Delta2inf)/(Psup-Pinf))*(P-Pinf);
+ result.nn11 = nn11inf + (nn11sup - nn11inf) / (Psup-Pinf) * (P-Pinf);
+ result.nn33 = nn33inf + (nn33sup - nn33inf) / (Psup-Pinf) * (P-Pinf);
index = k;
k=dataSize;
}
else if (data[k][1] == P)
- { //cerr << "elseif" << endl;
+ {
result.V = data[k][2];
result.F = data[k][3];
result.delta1 = data[k][4];
result.delta2 = data[k][5];
+ result.nn11= data[k][6];
+ result.nn33= data[k][7];
index = k;
k=dataSize;
=== modified file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp 2015-07-14 22:25:36 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp 2016-09-21 16:20:15 +0000
@@ -37,6 +37,8 @@
Real F; // adimentionnal capillary force for this meniscus : true force / ( 2 * pi * Rmax * superficial tension), (30) of Annexe1 of Scholtes2009d
Real delta1; // angle defined Fig 2.5 Scholtes2009d
Real delta2; // angle defined Fig 2.5 Scholtes2009d
+ Real nn11; // CapillaryPhys.nn11 / R2^2
+ Real nn33; // CapillaryPhys.nn33 / R2^2
int index1;
int index2;
@@ -94,6 +96,7 @@
((bool,binaryFusion,true,,"If true, capillary forces are set to zero as soon as, at least, 1 overlap (menisci fusion) is detected. Otherwise :yref:`fCap<CapillaryPhys.fCap>` = :yref:`fCap<CapillaryPhys.fCap>` / (:yref:`fusionNumber<CapillaryPhys.fusionNumber>` + 1 )"))
((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]_"))
((Real,surfaceTension,0.073,,"Value of considered surface tension")) // (0.073 N/m is water tension at 20 Celsius degrees)
+ ((string,suffCapFiles,"",,"Capillary files suffix: M(r=X)suffCapFiles"))
,,/*constructor*/
hertzInitialized = false;
hertzOn = false;
@@ -107,7 +110,7 @@
public:
Real D; // one cst D value in each TableauD (the one appearing last line of corresponding group D=cst in the capillary file)
std::vector<std::vector<Real> > data;
- MeniscusParameters Interpolate3(Real P, int& index);
+ MeniscusParameters Interpolate3(Real P, int& index); // does the interpolation on uc*
TableauD();
TableauD(std::ifstream& file);
~TableauD();
@@ -122,14 +125,14 @@
public:
Real R;
std::vector<TableauD> full_data; // members of full_data are the different TableauD, for different D.
- MeniscusParameters Interpolate2(Real D, Real P, int& index1, int& index2);
+ MeniscusParameters Interpolate2(Real D, Real P, int& index1, int& index2); // does the interpolation on d* (returning no meniscus when d* > the greatest D of the file)
std::ifstream& operator<< (std::ifstream& file);
Tableau();
Tableau(const char* filename);
~Tableau();
};
-class capillarylaw // the whole set of capillary files M(r=..)
+class capillarylaw // class for a whole set of capillary files M(r=..)
{
public:
capillarylaw();
