yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3403
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Thu 2014-02-06 19:48:43 +0100
message:
  add cell->info().trapCapP; fix pressure calculation for trapped phase.
modified:
  lib/triangulation/def_types.h
  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 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h	2014-01-24 18:37:37 +0000
+++ lib/triangulation/def_types.h	2014-02-06 18:48:43 +0000
@@ -226,11 +226,13 @@
 	Real capillaryCellVolume;//for calculating saturation
 	std::vector<double> poreRadius;//pore throat radius for drainage
 	double solidLine [4][4];//the length of intersecting line between sphere and facet. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_Lines in the facet.
+	double trapCapP;//for calculating the pressure of trapped phase, cell.Pressure_trapped = Pressure_air - trapCapP.
 	UnsatCellInfo (void)
 	{
 		poreRadius.resize(4, 0);
 		isWaterReservoir = false; isAirReservoir = false; capillaryCellVolume = 0;	  
-		for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidLine[k][l]=0;		
+		for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidLine[k][l]=0;
+		trapCapP = 0;
 	}
 };
 

=== modified file 'pkg/dem/UnsaturatedEngine.cpp'
--- pkg/dem/UnsaturatedEngine.cpp	2014-02-05 19:13:23 +0000
+++ pkg/dem/UnsaturatedEngine.cpp	2014-02-06 18:48:43 +0000
@@ -22,6 +22,7 @@
 
 #include "UnsaturatedEngine.hpp"
 #include <CGAL/Plane_3.h>
+#include <CGAL/Plane_3.h>
 
 CREATE_LOGGER ( UnsaturatedEngine );
 
@@ -149,10 +150,28 @@
 {
     updateReservoir(flow);
     RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
+    Finite_cells_iterator cell_end = tri.finite_cells_end();
+    for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+        if(cell->info().isAirReservoir == true)
+            invadeSingleCell2(cell,cell->info().p(),flow);
+    }
+    checkTrap(flow,bndCondValue[2]);
     Finite_cells_iterator _cell_end = tri.finite_cells_end();
     for ( Finite_cells_iterator _cell = tri.finite_cells_begin(); _cell != _cell_end; _cell++ ) {
-        if(_cell->info().isAirReservoir == true)
-            invadeSingleCell2(_cell,_cell->info().p(),flow);
+        if( (_cell->info().isWaterReservoir) || (_cell->info().isAirReservoir) ) continue;
+        _cell->info().p() = bndCondValue[2] - _cell->info().trapCapP;
+    }
+}
+
+template<class Solver>//check trapped phase, define trapCapP.
+void UnsaturatedEngine::checkTrap(Solver& flow, double pressure)
+{
+    RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
+    Finite_cells_iterator cell_end = tri.finite_cells_end();
+    for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+      if( (cell->info().isWaterReservoir) || (cell->info().isAirReservoir) ) continue;
+      if(cell->info().trapCapP!=0) continue;
+      cell->info().trapCapP=pressure;
     }
 }
 

=== modified file 'pkg/dem/UnsaturatedEngine.hpp'
--- pkg/dem/UnsaturatedEngine.hpp	2014-02-05 19:13:23 +0000
+++ pkg/dem/UnsaturatedEngine.hpp	2014-02-06 18:48:43 +0000
@@ -75,6 +75,7 @@
 		TPL void invadeSingleCell2(Cell_handle cell, double pressure, Solver& flow);
 		TPL void invade (Solver& flow );
 		TPL void invade2 (Solver& flow );
+		TPL void checkTrap(Solver& flow, double pressure);//check trapped phase, define trapCapP.		
 		TPL Real getMinEntryValue (Solver& flow );
 		TPL Real getMinEntryValue2 (Solver& flow);
 		TPL Real getSaturation(Solver& flow);
@@ -124,7 +125,7 @@
 		Real computePoreArea(Cellhandle cell, int j);
 		template<class Cellhandle>
 		Real computePorePerimeter(Cellhandle cell, int j);		
-		void saveVtk() {solver->saveVtk();}
+		void saveVtk() {bool initT= solver->noCache; solver->noCache=false; solver->saveVtk(); solver->noCache=initT;}
 		python::list getConstrictions() {
 			vector<Real> csd=solver->getConstrictions(); python::list pycsd;
 			for (unsigned int k=0;k<csd.size();k++) pycsd.append(csd[k]); return pycsd;}