yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #12129
[Branch ~yade-pkg/yade/git-trunk] Rev 3692: -set pore throat radius between two fictious cells negative in invadeBoundary=true mode.
------------------------------------------------------------
revno: 3692
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Tue 2015-06-30 16:20:30 +0200
message:
-set pore throat radius between two fictious cells negative in invadeBoundary=true mode.
modified:
pkg/pfv/TwoPhaseFlowEngine.cpp
pkg/pfv/TwoPhaseFlowEngine.hpp
pkg/pfv/UnsaturatedEngine.cpp
--
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/pfv/TwoPhaseFlowEngine.cpp'
--- pkg/pfv/TwoPhaseFlowEngine.cpp 2015-01-09 12:40:09 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.cpp 2015-06-30 14:20:30 +0000
@@ -437,6 +437,21 @@
}
vtkfile.end_data();
}
+void TwoPhaseFlowEngine::computePoreThroatRadiusTricky()
+{
+ computePoreThroatRadius();
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+ CellHandle neighbourCell;
+ for (FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++) {
+ for (int j=0; j<4; j++) {
+ neighbourCell = cell->neighbor(j);
+ if(cell->info().isFictious && neighbourCell->info().isFictious)
+ {cell->info().poreThroatRadius[j]=-1.0;
+ neighbourCell->info().poreThroatRadius[tri.mirror_index(cell, j)]= cell->info().poreThroatRadius[j];}
+ }
+ }
+}
void TwoPhaseFlowEngine::initialization()
{
@@ -444,7 +459,8 @@
setPositionsBuffer(true);//copy sphere positions in a buffer...
buildTriangulation(0.0,*solver);//create a triangulation and initialize pressure in the elements (connecting with W-reservoir), everything will be contained in "solver"
// initializeCellIndex();//initialize cell index
- computePoreThroatRadius();//save pore throat radius before drainage. Thomas, here you can also revert this to computePoreThroatCircleRadius().
+ if(isInvadeBoundary) {computePoreThroatRadius();}
+ else {computePoreThroatRadiusTricky();}//save pore throat radius before drainage. Thomas, here you can also revert this to computePoreThroatCircleRadius().
computePoreBodyRadius();//save pore body radius before imbibition
computePoreBodyVolume();//save capillary volume of all cells, for fast calculating saturation
computeSolidLine();//save cell->info().solidLine[j][y]
=== modified file 'pkg/pfv/TwoPhaseFlowEngine.hpp'
--- pkg/pfv/TwoPhaseFlowEngine.hpp 2015-05-29 16:53:54 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.hpp 2015-06-30 14:20:30 +0000
@@ -76,6 +76,8 @@
void computePoreCapillaryPressure(CellHandle cell);
void savePhaseVtk(const char* folder);
void computePoreThroatRadius();
+ void computePoreThroatRadiusTricky();//set the radius of pore throat between side pores negative.
+
double computeEffPoreThroatRadius(CellHandle cell, int j);
double computeEffPoreThroatRadiusFine(CellHandle cell, int j);
double bisection(CellHandle cell, int j, double a, double b);
@@ -122,6 +124,7 @@
((double,surfaceTension,0.0728,,"Water Surface Tension in contact with air at 20 Degrees Celsius is: 0.0728(N/m)"))
((bool,initialWetting,true,,"Initial wetting saturated (=true) or non-wetting saturated (=false)"))
((bool, isPhaseTrapped,true,,"If True, both phases can be entrapped by the other, which would correspond to snap-off. If false, both phases are always connected to their reservoirs, thus no snap-off."))
+ ((bool, isInvadeBoundary, true,,"Invasion side boundary condition. If True, pores of side boundary can be invaded; if False, the pore throats connecting side boundary are closed, those pores are excluded in saturation calculation."))
((bool, drainageFirst, true,,"If true, activate drainage first (initial saturated), then imbibition; if false, activate imbibition first (initial unsaturated), then drainage."))
((double,dtDynTPF,0.0,,"Parameter which stores the smallest time step, based on the residence time"))
=== modified file 'pkg/pfv/UnsaturatedEngine.cpp'
--- pkg/pfv/UnsaturatedEngine.cpp 2015-05-29 16:53:54 +0000
+++ pkg/pfv/UnsaturatedEngine.cpp 2015-06-30 14:20:30 +0000
@@ -20,7 +20,6 @@
// void initialization();
public :
- void computeTotalPoresVolume();
double computeCellInterfacialArea(CellHandle cell, int j, double rC);
// void computeSolidLine();
@@ -81,7 +80,6 @@
YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(UnsaturatedEngine,TwoPhaseFlowEngine,"Preliminary version engine of a drainage model for unsaturated soils. Note:Air reservoir is on the top; water reservoir is on the bottom.",
((bool, computeForceActivated, true,,"Activate capillary force computation. WARNING: turning off means capillary force is not computed at all, but the drainage can still work."))
- ((bool, isInvadeBoundary, true,,"Invasion side boundary condition. If True, pores of side boundary can be invaded; if False, the pore throats connecting side boundary are closed, those pores are excluded in saturation calculation."))
((int, windowsNo, 10,, "Number of genrated windows(or zoomed samples)."))
((bool, isDrainageActivated, true,, "Activates drainage."))
((bool, isImbibitionActivated, false,, "Activates imbibition."))
@@ -163,24 +161,6 @@
scene->forces.addForce ( V_it->info().id(), force); }}
*/}
-void UnsaturatedEngine::computeTotalPoresVolume()
-{
- initializeVolumes(*solver);
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- FiniteCellsIterator cellEnd = tri.finite_cells_end();
- totalCellVolume=0;
-
- if(isInvadeBoundary==true) {
- for (FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++) {
- if (cell->info().Pcondition) continue;//NOTE:reservoirs cells should not be included in totalCellVolume
- totalCellVolume = totalCellVolume + std::abs( cell->info().volume() );}}
- else {
- for (FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++) {
- if (cell->info().Pcondition) continue;//NOTE:reservoirs cells should not be included in totalCellVolume
- if (cell->info().isFictious) continue;
- totalCellVolume = totalCellVolume + std::abs( cell->info().volume() );}}
-}
-
void UnsaturatedEngine::updatePressure()
{
boundaryConditions(*solver);
@@ -194,8 +174,9 @@
if (isPhaseTrapped) {
if ( cell->info().isTrapW ) {cell->info().p()=bndCondValue[3]-cell->info().trapCapP;}
if ( cell->info().isTrapNW) {cell->info().p()=bndCondValue[2]+cell->info().trapCapP;}
- //here, in !isInvadeBoundary case, the pressure of bnd[0,1,4,5] cells is equal to the W-Res pressure, in order to calculate forces.
- if ( !cell->info().isWRes && !cell->info().isNWRes && !cell->info().isTrapW && !cell->info().isTrapNW ) {cell->info().p()=bndCondValue[2]; if (isInvadeBoundary) cerr<<"Something wrong in updatePressure.(isInvadeBoundary)";}
+ //check cell reservoir info.
+ if ( !cell->info().isWRes && !cell->info().isNWRes && !cell->info().isTrapW && !cell->info().isTrapNW ) {cerr<<"ERROR! NOT FIND Cell Info!";}
+// {cell->info().p()=bndCondValue[2]; if (isInvadeBoundary) cerr<<"Something wrong in updatePressure.(isInvadeBoundary)";}
}
}
}
@@ -215,7 +196,8 @@
CellHandle nCell = cell->neighbor(facet);
if (solver->T[solver->currentTes].Triangulation().is_infinite(nCell)) continue;
if (nCell->info().Pcondition) continue;//FIXME:defensive
- if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+// if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+ if (cell->info().poreThroatRadius[facet]<0) continue;
if ( (nCell->info().saturation==localSaturation) && (nCell->info().p() != localPressure) && ((nCell->info().isTrapNW)||(nCell->info().isTrapW)) ) {
nCell->info().p() = localPressure;
@@ -307,7 +289,7 @@
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
FiniteCellsIterator cellEnd = tri.finite_cells_end();
for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
- if( (cell->info().isFictious) && (!cell->info().Pcondition) && (!isInvadeBoundary) ) continue;
+// if( (cell->info().isFictious) && (!cell->info().Pcondition) && (!isInvadeBoundary) ) continue;
if( (cell->info().isWRes) || (cell->info().isNWRes) || (cell->info().isTrapW) || (cell->info().isTrapNW) ) continue;
cell->info().trapCapP=pressure;
if(cell->info().saturation==1.0) cell->info().isTrapW=true;
@@ -342,7 +324,7 @@
CellHandle nCell = cell->neighbor(facet);
if (solver->T[solver->currentTes].Triangulation().is_infinite(nCell)) continue;
if (nCell->info().Pcondition) continue;
- if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+// if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
if (nCell->info().saturation != 1.0) continue;
if (nCell->info().isWRes==true) continue;
nCell->info().isWRes = true;
@@ -359,7 +341,7 @@
CellHandle nCell = cell->neighbor(facet);
if (solver->T[solver->currentTes].Triangulation().is_infinite(nCell)) continue;
if (nCell->info().Pcondition) continue;
- if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+// if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
if (nCell->info().saturation != 0.0) continue;
if (nCell->info().isNWRes==true) continue;
nCell->info().isNWRes = true;
@@ -426,8 +408,8 @@
CellHandle nCell = cell->neighbor(facet);
if (tri.is_infinite(nCell)) continue;
if (nCell->info().Pcondition) continue;
- if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
- if ( nCell->info().isWRes == true ) {
+// if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+ if ( nCell->info().isWRes == true && cell->info().poreThroatRadius[facet]>0) {
double nCellP = std::max( (surfaceTension/cell->info().poreThroatRadius[facet]),(surfaceTension/nCell->info().poreBodyRadius) );
// double nCellP = surfaceTension/cell->info().poreThroatRadius[facet];
nextEntry = std::min(nextEntry,nCellP);}}}}
@@ -450,8 +432,8 @@
CellHandle nCell = cell->neighbor(facet);
if (tri.is_infinite(nCell)) continue;
if (nCell->info().Pcondition) continue;
- if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
- if ( nCell->info().isNWRes == true ) {
+// if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+ if ( nCell->info().isNWRes == true && cell->info().poreThroatRadius[facet]>0) {
double nCellP = std::min( (surfaceTension/nCell->info().poreBodyRadius), (surfaceTension/cell->info().poreThroatRadius[facet]));
nextEntry = std::max(nextEntry,nCellP);}}}}
