yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07420
[Branch ~yade-dev/yade/trunk] Rev 2806: - Fixed calculation of fluid velocity in cells
------------------------------------------------------------
revno: 2806
committer: Emanuele Catalano <catalano@xxxxxxxxxxx
branch nick: yade
timestamp: Wed 2011-04-06 18:37:32 +0200
message:
- Fixed calculation of fluid velocity in cells
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/def_types.h
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.hpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'lib/triangulation/FlowBoundingSphere.cpp'
--- lib/triangulation/FlowBoundingSphere.cpp 2011-04-06 12:17:47 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2011-04-06 16:37:32 +0000
@@ -299,7 +299,7 @@
return Tes;
}
-void FlowBoundingSphere::Average_Cell_Velocity()
+void FlowBoundingSphere::Average_Relative_Cell_Velocity()
{
RTriangulation& Tri = T[currentTes].Triangulation();
Point pos_av_facet;
@@ -323,8 +323,7 @@
pos_av_facet = (Point) cell->info() + ( branch*Surfk ) *Surfk;
// pos_av_facet=CGAL::ORIGIN + ((cell->vertex(facetVertices[i][0])->point() - CGAL::ORIGIN) + (cell->vertex(facetVertices[i][1])->point() - CGAL::ORIGIN) + (cell->vertex(facetVertices[i][2])->point() - CGAL::ORIGIN))*0.3333333333;
facet_flow_rate = (cell->info().k_norm())[i] * (cell->info().p() - cell->neighbor (i)->info().p());
- for (int y=0;y<3;y++) volume_facet_translation += (cell->info().facetVelocity())[i][y]*cell->info().facetSurfaces[i][y];
- cell->info().av_vel() = cell->info().av_vel() + (facet_flow_rate - volume_facet_translation) * ( pos_av_facet-CGAL::ORIGIN );
+ cell->info().av_vel() = cell->info().av_vel() + (facet_flow_rate) * ( pos_av_facet-CGAL::ORIGIN );
}}
if (cell->info().volume()){ tVel+=cell->info().av_vel()[1]; tVol+=cell->info().volume();}
cell->info().av_vel() = cell->info().av_vel() /cell->info().volume();
@@ -347,7 +346,7 @@
void FlowBoundingSphere::Average_Fluid_Velocity()
{
- Average_Cell_Velocity();
+ Average_Relative_Cell_Velocity();
RTriangulation& Tri = T[currentTes].Triangulation();
int num_vertex = 0;
Finite_vertices_iterator vertices_end = Tri.finite_vertices_end();
@@ -393,7 +392,7 @@
vector<Real> FlowBoundingSphere::Average_Fluid_Velocity_On_Sphere(int Id_sph)
{
- Average_Cell_Velocity();
+ Average_Relative_Cell_Velocity();
RTriangulation& Tri = T[currentTes].Triangulation();
Real Volumes; CGT::Vecteur VelocityVolumes;
@@ -1526,7 +1525,7 @@
void FlowBoundingSphere::ComsolField()
{
//Compute av. velocity first, because in the following permeabilities will be overwritten with "junk" (in fact velocities from comsol)
- Average_Cell_Velocity();
+ Average_Relative_Cell_Velocity();
RTriangulation& Tri = T[currentTes].Triangulation();
Cell_handle c;
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2011-04-05 16:58:41 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2011-04-06 16:37:32 +0000
@@ -110,7 +110,7 @@
void ComsolField();
void Interpolate ( Tesselation& Tes, Tesselation& NewTes );
- void Average_Cell_Velocity();
+ void Average_Relative_Cell_Velocity();
void Average_Fluid_Velocity();
void ApplySinusoidalPressure(RTriangulation& Tri, double Amplitude, double Average_Pressure, double load_intervals);
bool isOnSolid (double X, double Y, double Z);
=== modified file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h 2011-04-06 12:17:47 +0000
+++ lib/triangulation/def_types.h 2011-04-06 16:37:32 +0000
@@ -54,13 +54,12 @@
int fict;
Real VolumeVariation;
double pression; //stockage d'une valeur de pression pour chaque cellule
- Vecteur Average_Cell_Velocity; //average velocity defined for a single cell as 1/Volume * SUM_ON_FACETS(x_average_facet*average_facet_flow_rate)
+ Vecteur Average_Cell_Velocity; //average relative (fluid - facet translation) velocity defined for a single cell as 1/Volume * SUM_ON_FACETS(x_average_facet*average_facet_flow_rate)
// Surface vectors of facets, pointing from outside toward inside the cell
std::vector<Vecteur> facetSurfaces;
//Ratio between fluid surface and facet surface
std::vector<Real> facetFluidSurfacesRatio;
- std::vector<Vecteur> facetVelocities;
// Reflects the geometrical property of the cell, so that the force by cell fluid on grain "i" is pressure*unitForceVectors[i]
std::vector<Vecteur> unitForceVectors;
// Store the area of triangle-sphere intersections for each facet (used in forces definition)
@@ -79,7 +78,6 @@
facetSurfaces.resize(4);
facetFluidSurfacesRatio.resize(4);
facetSphereCrossSections.resize(4);
- facetVelocities.resize(4);
unitForceVectors.resize(4);
for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidSurfaces[k][l]=0;
RayHydr.resize(4, 0);
@@ -119,7 +117,6 @@
inline std::vector<double>& k_norm (void) {return module_permeability;}
inline std::vector< Vecteur >& facetSurf (void) {return facetSurfaces;}
- inline std::vector< Vecteur >& facetVelocity (void) {return facetVelocities;}
inline std::vector<Vecteur>& force (void) {return cell_force;}
inline std::vector<double>& Rh (void) {return RayHydr;}
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2011-04-05 16:58:41 +0000
+++ pkg/dem/FlowEngine.cpp 2011-04-06 16:37:32 +0000
@@ -21,12 +21,6 @@
CREATE_LOGGER (FlowEngine);
-BOOST_PYTHON_MODULE(_vtkFile){
- python::scope().attr("__doc__")="SaveVTKFile";
- YADE_SET_DOCSTRING_OPTS;
- python::class_<FlowEngine>("FlowEngine","Create file for parallel visualization" )
- .def("saveVtk",&FlowEngine::saveVtk,"Save VTK File, each dd iterations");}
-
FlowEngine::~FlowEngine()
{
}
@@ -210,6 +204,7 @@
flow->DEBUG_OUT = Debug;
flow->useSolver = useSolver;
flow->VISCOSITY = viscosity;
+ flow->areaR2Permeability=areaR2Permeability;
flow->T[flow->currentTes].Clear();
flow->T[flow->currentTes].max_id=-1;
@@ -392,6 +387,64 @@
if (Debug) cout << "Volumes initialised." << endl;
}
+void FlowEngine::Average_real_cell_velocity()
+{
+ flow->Average_Relative_Cell_Velocity();
+ Vector3r Vel (0,0,0);
+ //AVERAGE CELL VELOCITY
+ CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
+ for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
+ switch ( cell->info().fictious()) {
+ case ( 3 ):
+ for ( int g=0;g<4;g++ )
+ {
+ if ( !cell->vertex ( g )->info().isFictious ) {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ for (int i=0;i<3;i++) Vel[i] = Vel[i] + sph->state->vel[i]/4;}
+ }
+ break;
+ case ( 2 ):
+ for ( int g=0;g<4;g++ )
+ {
+ if ( !cell->vertex ( g )->info().isFictious ) {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ for (int i=0;i<3;i++) Vel[i] = Vel[i] + sph->state->vel[i]/4;}
+ }
+ break;
+ case ( 1 ):
+ for ( int g=0;g<4;g++ )
+ {
+ if ( !cell->vertex ( g )->info().isFictious ) {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ for (int i=0;i<3;i++) Vel[i] = Vel[i] + sph->state->vel[i]/4;}
+ }
+ break;
+ case ( 0 ) :
+ for ( int g=0;g<4;g++ )
+ {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ for (int i=0;i<3;i++) Vel[i] = Vel[i] + sph->state->vel[i]/4;}
+ }
+ break;
+
+
+ CGT::RTriangulation& Tri = flow->T[flow->currentTes].Triangulation();
+ CGT::Point pos_av_facet;
+ double volume_facet_translation = 0;
+ CGT::Vecteur Vel_av (Vel[0], Vel[1], Vel[2]);
+ for ( int i=0; i<4; i++ ) {
+ volume_facet_translation = 0;
+ if (!Tri.is_infinite(cell->neighbor(i))){
+ CGT::Vecteur Surfk = cell->info()-cell->neighbor(i)->info();
+ Real area = sqrt ( Surfk.squared_length() );
+ Surfk = Surfk/area;
+ CGT::Vecteur branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
+ pos_av_facet = (CGT::Point) cell->info() + ( branch*Surfk ) *Surfk;
+ volume_facet_translation += Vel_av*cell->info().facetSurfaces[i];
+ cell->info().av_vel() = cell->info().av_vel() - volume_facet_translation/cell->info().volume() * ( pos_av_facet-CGAL::ORIGIN );}}
+ }
+}
+
void FlowEngine::UpdateVolumes ()
{
if (Debug) cout << "Updating volumes.............." << endl;
@@ -429,10 +482,6 @@
int w=0;
Real Wall_coordinate=0;
-
- Vector3r Vel[3];
- int id = 0;
- double Vel_x=0, Vel_y=0, Vel_z=0;
for ( int y=0;y<4;y++ )
{
@@ -441,7 +490,6 @@
const shared_ptr<Body>& sph = Body::byId
( cell->vertex ( y )->info().id(), scene );
for ( int g=0;g<3;g++ ) V[w][g]=sph->state->pos[g];
- Vel[w] = sph->state->pos;
w++;
}
else
@@ -450,23 +498,8 @@
const shared_ptr<Body>& wll = Body::byId ( b , scene );
if (!flow->boundaries[b].useMaxMin) Wall_coordinate = wll->state->pos[flow->boundaries[b].coordinate]+(flow->boundaries[b].normal[flow->boundaries[b].coordinate])*wall_thickness/2;
else Wall_coordinate = flow->boundaries[b].p[flow->boundaries[b].coordinate];
- id = y;
}
}
-
- for ( int y=0;y<4;y++ ){
- for ( int j=0;j<3;j++ ){
- if (cell->vertex(j)->info().isFictious){
- Vel_x += 0.33333333333*(Vel[facetVertices[y][j]][0]);
- Vel_y += 0.33333333333*(Vel[facetVertices[y][j]][1]);
- Vel_z += 0.33333333333*(Vel[facetVertices[y][j]][2]);}
- else for (int j2=0;j2<3;j2++){
- if (!cell->vertex(j2)->info().isFictious){
- Vel_x += 0.5*(Vel[facetVertices[y][j]][0]);
- Vel_y += 0.5*(Vel[facetVertices[y][j]][1]);
- Vel_z += 0.5*(Vel[facetVertices[y][j]][2]);}}}
- CGT::Vecteur Vel_facet ( Vel_x, Vel_y, Vel_z );
- (cell->info().facetVelocity())[y] = Vel_facet*1;}
double v1[3], v2[3];
@@ -490,10 +523,6 @@
int j=0;
bool first_sph=true;
- Vector3r Vel[2];
- vector<int> id;
- id.resize(2);
-
for ( int g=0;g<4;g++ )
{
if ( cell->vertex ( g )->info().isFictious )
@@ -509,24 +538,15 @@
const shared_ptr<Body>& sph1 = Body::byId
( cell->vertex ( g )->info().id(), scene );
for ( int k=0;k<3;k++ ) { A[k]=AS[k]=AT[k]=sph1->state->pos[k]; first_sph=false;}
- Vel[0] = sph1->state->vel; id[0]=g;
}
else
{
const shared_ptr<Body>& sph2 = Body::byId
( cell->vertex ( g )->info().id(), scene );
for ( int k=0;k<3;k++ ) { B[k]=BS[k]=BT[k]=sph2->state->pos[k]; }
- Vel[1] = sph2->state->vel; id[1]=g;
}
}
- for (int y=0;y<4;y++){
- if ( y == id[0] ) {CGT::Vecteur Vel_facet (Vel[0][0], Vel[0][1], Vel[0][2]); (cell->info().facetVelocity())[y] = Vel_facet*1;}
- else if ( y == id[1] ) {CGT::Vecteur Vel_facet (Vel[1][0], Vel[1][1], Vel[1][2]); (cell->info().facetVelocity())[y] = Vel_facet*1;}
- else { CGT::Vecteur Vel_facet (0.5*(Vel[0][0]+Vel[0][1]+Vel[0][2]),0.5*(Vel[1][0]+Vel[1][1]+Vel[1][2]),0.5*(Vel[2][0]+Vel[2][1]+Vel[2][2]));
- (cell->info().facetVelocity())[y] = Vel_facet*1;}}
-
-
AS[flow->boundaries[b[0]].coordinate]=BS[flow->boundaries[b[0]].coordinate] = Wall_coordinate[0];
AT[flow->boundaries[b[1]].coordinate]=BT[flow->boundaries[b[1]].coordinate] = Wall_coordinate[1];
@@ -550,8 +570,6 @@
Real Wall_coordinate[3];
int j=0;
- Vector3r Vel;
-
for ( int g=0;g<4;g++ )
{
if ( cell->vertex ( g )->info().isFictious )
@@ -566,13 +584,9 @@
{
const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
for ( int k=0;k<3;k++ ) { A[k]=AS[k]=AT[k]=AW[k]=sph->state->pos[k];}
- Vel = sph->state->vel;
}
}
- CGT::Vecteur Vel_facet (Vel[0], Vel[1], Vel[2]);
- for (int y=0;y<4;y++) (cell->info().facetVelocity())[y] = Vel_facet*1;
-
AS[flow->boundaries[b[0]].coordinate]= AT[flow->boundaries[b[0]].coordinate]= AW[flow->boundaries[b[0]].coordinate]= Wall_coordinate[0];
AT[flow->boundaries[b[1]].coordinate]= Wall_coordinate[1];
AW[flow->boundaries[b[2]].coordinate]= Wall_coordinate[2];
@@ -590,24 +604,6 @@
Real FlowEngine::Volume_cell ( CGT::Cell_handle cell)
{
Vector3r A[4];
- Vector3r Vel[4];
-
- double Vel_x=0, Vel_y=0, Vel_z=0;
-
- for ( int y=0;y<4;y++ )
- {
- const shared_ptr<Body>& sph = Body::byId(cell->vertex ( y )->info().id(), scene);
- A[y]=sph->state->pos;
- Vel[y]=sph->state->vel;
- }
-
- for ( int y=0;y<4;y++ ){
- for ( int j=0;j<3;j++ ){
- Vel_x += 0.33333333333*(Vel[facetVertices[y][j]][0]);
- Vel_y += 0.33333333333*(Vel[facetVertices[y][j]][1]);
- Vel_z += 0.33333333333*(Vel[facetVertices[y][j]][2]);}
- CGT::Vecteur Vel_facet ( Vel_x, Vel_y, Vel_z );
- (cell->info().facetVelocity())[y] = Vel_facet*1; }
CGT::Point p1 ( ( A[0] ) [0], ( A[0] ) [1], ( A[0] ) [2] );
CGT::Point p2 ( ( A[1] ) [0], ( A[1] ) [1], ( A[1] ) [2] );
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2011-04-05 16:58:41 +0000
+++ pkg/dem/FlowEngine.hpp 2011-04-06 16:37:32 +0000
@@ -51,6 +51,7 @@
void BoundaryConditions();
void imposePressure(Vector3r pos, Real p);
void clearImposedPressure();
+ void Average_real_cell_velocity();
Real getFlux(int cond);
void saveVtk() {flow->save_vtk_file();}
vector<Real> AvFlVelOnSph(int id_sph) {flow->Average_Fluid_Velocity_On_Sphere(id_sph);}
@@ -100,6 +101,7 @@
((int, intervals, 30,, "Number of layers for computation average fluid pressure profiles to build consolidation curves"))
((int, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso"))
((bool, liquefaction, false,,"Measure fluid pressure along the heigth of the sample"))
+ ((double, V_d, 0,,"darcy velocity of fluid in sample"))
// ((double, bottom_seabed_pressure,0,,"Fluid pressure measured at the bottom of the seabed on the symmetry axe"))
((bool, Flow_imposed_TOP_Boundary, true,, "if false involve pressure imposed condition"))
((bool, Flow_imposed_BOTTOM_Boundary, true,, "if false involve pressure imposed condition"))
@@ -120,6 +122,7 @@
((bool, LEFT_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
((bool, FRONT_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
((bool, BACK_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
+ ((bool, areaR2Permeability, 1,,"Use corrected formula for permeabilities calculation in flowboundingsphere (areaR2permeability variable)"))
,,
timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
,
