yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11411
[Branch ~yade-pkg/yade/git-trunk] Rev 3407: add getSaturation2() for mode 2.
------------------------------------------------------------
revno: 3407
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Thu 2014-02-20 18:11:45 +0100
message:
add getSaturation2() for mode 2.
modified:
pkg/dem/UnsaturatedEngine.cpp
pkg/dem/UnsaturatedEngine.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/UnsaturatedEngine.cpp'
--- pkg/dem/UnsaturatedEngine.cpp 2014-02-19 11:07:57 +0000
+++ pkg/dem/UnsaturatedEngine.cpp 2014-02-20 17:11:45 +0000
@@ -56,25 +56,32 @@
}
void UnsaturatedEngine::action()
-{
- if ( !isActivated ) return;
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- if ( (tri.number_of_vertices()==0) || (solver->updateTriangulation) ) {
- cout<< "triangulation is empty: building a new one" << endl;
- scene = Omega::instance().getScene().get();//here define the pointer to Yade's scene
- setPositionsBuffer(true);//copy sphere positions in a buffer...
- Build_Triangulation(P_zero,solver);//create a triangulation and initialize pressure in the elements, everything will be contained in "solver"
- initializeCellIndex(solver);//initialize cell index
- initializePoreRadius(solver);//save all pore radii before invade
- updateVolumeCapillaryCell(solver);//save capillary volume of all cells, for calculating saturation
- computeSolidLine(solver);//save cell->info().solidLine[j][y]
- solver->noCache = true;
- }
- ///compute invade
- /*if ( solver->pressureChanged )*/ solver->invade();
- ///compute force
- solver->computeFacetPoreForcesWithCache();
- scene->forces.addForce ( V_it->info().id(), force);
+{/*
+ if ( !isActivated ) return;
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ if ( (tri.number_of_vertices()==0) || (updateTriangulation) ) {
+ cout<< "triangulation is empty: building a new one" << endl;
+ scene = Omega::instance().getScene().get();//here define the pointer to Yade's scene
+ setPositionsBuffer(true);//copy sphere positions in a buffer...
+ Build_Triangulation(P_zero,solver);//create a triangulation and initialize pressure in the elements, everything will be contained in "solver"
+ initializeCellIndex(solver);//initialize cell index
+ initializePoreRadius(solver);//save all pore radii before invade
+ updateVolumeCapillaryCell(solver);//save capillary volume of all cells, for calculating saturation
+ computeSolidLine(solver);//save cell->info().solidLine[j][y]
+ solver->noCache = true;
+ }
+ ///compute invade
+ if (pressureForce) {
+ invade(solver);
+ }
+ ///compute force
+ computeFacetPoreForcesWithCache(solver);
+ Vector3r force;
+ Finite_vertices_iterator vertices_end = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
+ for ( Finite_vertices_iterator V_it = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); V_it != vertices_end; V_it++ ) {
+ force = pressureForce ? Vector3r ( V_it->info().forces[0],V_it->info().forces[1],V_it->info().forces[2] ): Vector3r(0,0,0);
+ scene->forces.addForce ( V_it->info().id(), force);
+ }*/
}
///invade mode 1. update phase reservoir before invasion. Consider no viscous effects, and invade gradually.
@@ -146,7 +153,8 @@
}
}
if (nextEntry==1e50) {
- cout << "please set initial air pressure for the cell !" << endl;
+ cout << "End drainage !" << endl;
+ return nextEntry=0;
}
else {
return nextEntry;
@@ -263,6 +271,27 @@
}
}
+template<class Solver>
+Real UnsaturatedEngine::getSaturation (Solver& flow )
+{
+ updateReservoir(flow);
+ RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
+ Real capillary_volume = 0.0; //total capillary volume
+ Real air_volume = 0.0; //air volume
+ Finite_cells_iterator cell_end = tri.finite_cells_end();
+ for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+ if (tri.is_infinite(cell)) continue;
+ if (cell->info().Pcondition) continue;//when calculating saturation, exclude boundary cells?(chao)
+// if (cell.has_vertex() )
+ capillary_volume = capillary_volume + cell->info().capillaryCellVolume;
+ if (cell->info().isAirReservoir==true) {
+ air_volume = air_volume + cell->info().capillaryCellVolume;
+ }
+ }/*cerr<<"air_volume:"<<air_volume<<" capillary_volume:"<<capillary_volume<<endl;*/
+ Real saturation = 1 - air_volume/capillary_volume;
+ return saturation;
+}
+
///invade mode 2. Consider no trapped phase.
template<class Solver>
void UnsaturatedEngine::invadeSingleCell2(Cell_handle cell, double pressure, Solver& flow)
@@ -333,6 +362,26 @@
}
}
+template<class Solver>
+Real UnsaturatedEngine::getSaturation2(Solver& flow )
+{
+ RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
+ Real capillary_volume = 0.0;
+ Real water_volume = 0.0;
+ Finite_cells_iterator cell_end = tri.finite_cells_end();
+ for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+ if (tri.is_infinite(cell)) continue;
+ if (cell->info().Pcondition) continue;//when calculating saturation, exclude boundary cells?(chao)
+// if (cell.has_vertex() )
+ capillary_volume = capillary_volume + cell->info().capillaryCellVolume;
+ if (cell->info().p()==0) {
+ water_volume = water_volume + cell->info().capillaryCellVolume;
+ }
+ }
+ Real saturation = water_volume/capillary_volume;
+ return saturation;
+}
+
template<class Cellhandle>
double UnsaturatedEngine::computeEffPoreRadius(Cellhandle cell, int j)
{
@@ -812,27 +861,6 @@
}
template<class Solver>
-Real UnsaturatedEngine::getSaturation (Solver& flow )
-{
- updateReservoir(flow);
- RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
- Real capillary_volume = 0.0; //total capillary volume
- Real air_volume = 0.0; //air volume
- Finite_cells_iterator cell_end = tri.finite_cells_end();
- for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
- if (tri.is_infinite(cell)) continue;
- if (cell->info().Pcondition) continue;//when calculating saturation, exclude boundary cells?(chao)
-// if (cell.has_vertex() )
- capillary_volume = capillary_volume + cell->info().capillaryCellVolume;
- if (cell->info().isAirReservoir==true) {
- air_volume = air_volume + cell->info().capillaryCellVolume;
- }
- }/*cerr<<"air_volume:"<<air_volume<<" capillary_volume:"<<capillary_volume<<endl;*/
- Real saturation = 1 - air_volume/capillary_volume;
- return saturation;
-}
-
-template<class Solver>
void UnsaturatedEngine::saveListNodes(Solver& flow)
{
RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
=== modified file 'pkg/dem/UnsaturatedEngine.hpp'
--- pkg/dem/UnsaturatedEngine.hpp 2014-02-19 11:07:57 +0000
+++ pkg/dem/UnsaturatedEngine.hpp 2014-02-20 17:11:45 +0000
@@ -78,7 +78,9 @@
TPL void checkTrap(Solver& flow, double pressure);//check trapped phase, define trapCapP.
TPL Real getMinEntryValue (Solver& flow );
TPL Real getMinEntryValue2 (Solver& flow);
+ TPL void updatePressure2(Solver& flow);
TPL Real getSaturation(Solver& flow);
+ TPL Real getSaturation2(Solver& flow);
TPL void saveListNodes(Solver& flow);
TPL void saveListConnection(Solver& flow);
@@ -141,6 +143,7 @@
Real _getMinEntryValue() {return getMinEntryValue(solver);}
Real _getMinEntryValue2() {return getMinEntryValue2(solver);}
Real _getSaturation () {return getSaturation(solver);}
+ Real _getSaturation2() {return getSaturation2(solver);}
void _saveListNodes() {saveListNodes(solver);}
void _saveListConnection() {saveListConnection(solver);}
void _saveLatticeNodeX(double x) {saveLatticeNodeX(solver,x);}
@@ -164,6 +167,7 @@
((bool, Debug, false,,"Activate debug messages"))
((double, wall_thickness,0.001,,"Walls thickness"))
((double,P_zero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model."))
+ ((bool, updateTriangulation, 0,,"If true the medium is retriangulated."))
((double,gasPressure,0,,"Invasion pressure"))
((double,surfaceTension,0.0728,,"Water Surface Tension in contact with air at 20 Degrees Celsius is: 0.0728(N/m)"))
((double, porosity, 0,,"Porosity computed at each retriangulation"))
@@ -202,6 +206,7 @@
.def("testFunction",&UnsaturatedEngine::testFunction,"The playground for Chao's experiments.")
.def("buildTriangulation",&UnsaturatedEngine::_buildTriangulation,"Triangulate spheres of the current scene.")
.def("getSaturation",&UnsaturatedEngine::_getSaturation,"get saturation")
+ .def("getSaturation2",&UnsaturatedEngine::_getSaturation2,"get saturation for mode 2")
.def("getMinEntryValue",&UnsaturatedEngine::_getMinEntryValue,"get the minimum air entry pressure for the next invade step")
.def("getMinEntryValue2",&UnsaturatedEngine::_getMinEntryValue2,"get the minimum air entry pressure for the next invade step(version2)")
.def("saveListNodes",&UnsaturatedEngine::_saveListNodes,"Save the list of nodes.")
