------------------------------------------------------------ revno: 3934 committer: jduriez <jerome.dur...@ucalgary.ca> 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();
_______________________________________________ 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