yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07413
[Branch ~yade-dev/yade/trunk] Rev 2804: - re-fix fluid area definition
------------------------------------------------------------
revno: 2804
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Tue 2011-04-05 18:58:41 +0200
message:
- re-fix fluid area definition
- some cleaning here and there (Donia+Bruno)
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
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-01 09:59:17 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2011-04-05 16:58:41 +0000
@@ -67,7 +67,6 @@
FlowBoundingSphere::FlowBoundingSphere()
{
- id_Sphere=0;
x_min = 1000.0, x_max = -10000.0, y_min = 1000.0, y_max = -10000.0, z_min = 1000.0, z_max = -10000.0;
currentTes = 0;
nOfSpheres = 0;
@@ -79,10 +78,9 @@
for (int i=0;i<6;i++) boundsIds[i] = 0;
minPermLength=-1;
SLIP_ON_LATERALS = false;//no-slip/symmetry conditions on lateral boundaries
- TOLERANCE = 1e-06;
+ TOLERANCE = 1e-07;
RELAX = 1.9;
ks=0;
- V_darcy_Donia=0;
distance_correction = true;
meanK_LIMIT = true;
meanK_STAT = false; K_opt_factor=0;
@@ -92,7 +90,7 @@
RAVERAGE = false; /** use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
OUTPUT_BOUDARIES_RADII = false;
RAVERAGE = false; /** if true use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
- areaR2Permeability=true;
+ areaR2Permeability=false;
}
void FlowBoundingSphere::ResetNetwork() {noCache=true;}
@@ -686,9 +684,7 @@
Cell_handle neighbour_cell;
double k=0, distance = 0, radius = 0, viscosity = VISCOSITY;
- double m3=0, m1=0, m2=0, d=0, h=0;
int surfneg=0;
- Real S0=0;
int NEG=0, POS=0, pass=0;
bool ref = Tri.finite_cells_begin()->info().isvisited;
@@ -713,7 +709,6 @@
Sphere& v0 = W[0]->point();
Sphere& v1 = W[1]->point();
Sphere& v2 = W[2]->point();
- Sphere& t0=v0, t1=v1, t2=v2;
#ifdef USE_FAST_MATH
//FIXME : code not compiling,, do the same as in "else"
assert((W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point())>=0 && (W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point())<=1);
@@ -736,7 +731,7 @@
if (radius==0) {
cout << "INS-INS PROBLEM!!!!!!!" << endl;
}
-
+ Real h,d;
if (distance!=0) {
if (minPermLength>0 && distance_correction) distance=max(minPermLength,distance);
Real fluidArea=0;
@@ -744,54 +739,42 @@
Real area = sqrt(Surfk.squared_length());
const Vecteur& crossSections = cell->info().facetSphereCrossSections[j];
if (areaR2Permeability){
- //if (cell->info().fictious()==0 && neighbour_cell->info().fictious()==0){
-
- m1=sqrt((cross_product((v0-v1),v2-v1)).squared_length()/(v2-v1).squared_length());
- m2=sqrt((cross_product((v0-v1),v2-v0)).squared_length()/(v2-v0).squared_length());
- m3=sqrt((cross_product((v0-v2),v1-v0)).squared_length()/(v1-v0).squared_length());
-
- if (m1<sqrt(v0.weight()))
- {
- d=2*sqrt((v0.weight()-m1*m1));
- h=sqrt(v0.weight())-m1;
- S0=0.25*M_PI*(d*d+4*h*h);}
-
- if (m2<sqrt(v1.weight()))
- {
- d=2*sqrt((v1.weight()-m2*m2));
- h=sqrt(v1.weight())-m2;
- S0=0.25*M_PI*(d*d+4*h*h);}
-//
- if (m3<sqrt(v2.weight()))
- {
- d=2*sqrt((v2.weight()-m3*m3));
- h=sqrt(v2.weight())-m3;
- S0=0.25*M_PI*(d*d+4*h*h);}
-
-// if (S0>0) cout<<"S0= "<<S0<<endl;
-
- // }
+ Real m1=sqrt((cross_product((v0-v1),v2-v1)).squared_length()/(v2-v1).squared_length());
+ Real S0=0;
+ if (m1<sqrt(v0.weight())) {
+ d=2*sqrt((v0.weight()-m1*m1));
+ h=sqrt(v0.weight())-m1;
+ S0=0.25*M_PI*(d*d+4*h*h);}
+ else {
+ Real m2=sqrt((cross_product((v0-v1),v2-v0)).squared_length()/(v2-v0).squared_length());
+ if (m2<sqrt(v1.weight())) {
+ d=2*sqrt((v1.weight()-m2*m2));
+ h=sqrt(v1.weight())-m2;
+ S0=0.25*M_PI*(d*d+4*h*h);}
+ else {
+ Real m3=sqrt((cross_product((v0-v2),v1-v0)).squared_length()/(v1-v0).squared_length());
+ if (m3<sqrt(v2.weight())) {
+ d=2*sqrt((v2.weight()-m3*m3));
+ h=sqrt(v2.weight())-m3;
+ S0=0.25*M_PI*(d*d+4*h*h);}
+ }
+ }
fluidArea=area-crossSections[0]-crossSections[1]-crossSections[2]+S0;
k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);}
- else k = (M_PI * pow(radius,4)) / (8*viscosity*distance);
-
+ else k = (M_PI * pow(radius,4)) / (8*viscosity*distance);
+
if (k<0) {surfneg+=1;
cout<<"__ k<0 __"<<k<<" "<<area<<" "<<crossSections[0]<<" "<<crossSections[1]<<" "<<crossSections[2] <<endl;}
} else k = infiniteK;//Will be corrected in the next loop
(cell->info().k_norm())[j]= k*k_factor;
-// (cell->info().facetSurf())[j]= k*n;
(neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]= k*k_factor;
meanDistance += distance;
meanRadius += radius;
meanK += (cell->info().k_norm())[j];
-
-// if (!meanK_LIMIT) kFile << ( cell->info().k_norm() )[j] << endl;
-// (neighbour_cell->info().facetSurf())[Tri.mirror_index(cell, j)]= (-k) *n;
}
- // else if ( Tri.is_infinite ( neighbour_cell )) k = 0;//connection to an infinite cells
}
cell->info().isvisited = !ref;
}
@@ -1115,6 +1098,7 @@
// cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;}
// // save_vtk_file ( Tri );
// }
+// if (DEBUG_OUT) {cout << "pmax " << p_max << "; pmoy : " << p_moy << endl; cout << "iteration " << j <<"; erreur : " << dp_max/p_max << " tolerance " << tolerance<<endl;}
#ifdef GS_OPEN_MP
} while (j<1500);
#else
@@ -1197,13 +1181,8 @@
double viscosity = VISCOSITY;
double gravity = 9.80665;
double Vdarcy = Q1/Section;
- V_darcy_Donia=Vdarcy;
-// double GradP = abs(P_Inf-P_Sup) /DeltaY;
double DeltaP = abs(P_Inf-P_Sup);
-// double GradH = GradP/ (density*gravity);
double DeltaH = DeltaP/ (density*gravity);
-// double Ks= (Vdarcy) /GradH;
-// double k = Ks*viscosity/ (density*gravity);
double k = viscosity*Vdarcy*DeltaY / DeltaP; /**m²**/
double Ks = k*(density*gravity)/viscosity; /**m/s**/
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2011-04-01 09:59:17 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2011-04-05 16:58:41 +0000
@@ -45,7 +45,6 @@
bool OUTPUT_BOUDARIES_RADII;
bool noCache;//flag for checking if cached values cell->unitForceVectors have been defined
vector<pair<Point,Real> > imposedP;
- double V_darcy_Donia;
void initNewTri () {noCache=true; /*isLinearSystemSet=false; areCellsOrdered=false;*/}//set flags after retriangulation
bool computeAllCells;//exececute computeHydraulicRadius for all facets and all spheres (double cpu time but needed for now in order to define crossSections correctly)
@@ -54,9 +53,6 @@
bool RAVERAGE;
int walls_id[6];
-
- int id_Sphere;
-
void mplot ( char *filename);
void Localize();
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2011-03-23 09:38:16 +0000
+++ pkg/dem/FlowEngine.cpp 2011-04-05 16:58:41 +0000
@@ -93,8 +93,8 @@
{
id = V_it->info().id();
for (int y=0;y<3;y++) f[y] = (V_it->info().forces)[y];
- if (id==id_sphere)
- id_force = f;
+// if (id==id_sphere)
+// id_force = f;
// 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++ ){
// for (int j=0;j<4;j++){
@@ -210,8 +210,6 @@
flow->DEBUG_OUT = Debug;
flow->useSolver = useSolver;
flow->VISCOSITY = viscosity;
-
- flow->id_Sphere=id_sphere;
flow->T[flow->currentTes].Clear();
flow->T[flow->currentTes].max_id=-1;
@@ -233,28 +231,30 @@
if (first)
{
+ flow->TOLERANCE=Tolerance;
+ flow->RELAX=Relax;
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++ ){cell->info().dv() = 0; cell->info().p() = 0;}
- if (compute_K) {flow->TOLERANCE=1e-06; K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
+ if (compute_K) {K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
BoundaryConditions();
flow->Initialize_pressures( P_zero );
if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Amplitude, Sinus_Average, 30);
- flow->TOLERANCE=Tolerance;
- flow->RELAX=Relax;
+
+
}
else
{
+ flow->TOLERANCE=Tolerance;
+ flow->RELAX=Relax;
if (Debug && compute_K) cout << "---------UPDATE PERMEABILITY VALUE--------------" << endl;
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++ ){cell->info().dv() = 0; cell->info().p() = 0;}
- if (compute_K) {flow->TOLERANCE=1e-06; K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
+ if (compute_K) {K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
BoundaryConditions();
flow->Initialize_pressures(P_zero);// FIXME : why, if we are going to interpolate after that?
- flow->TOLERANCE=Tolerance;//So it can be changed at run time
flow->Interpolate (flow->T[!flow->currentTes], flow->T[flow->currentTes]);
if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Amplitude, Sinus_Average, 30);
- flow->TOLERANCE=Tolerance;
- flow->RELAX=Relax;
+
}
Initialize_volumes();
}
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2011-03-23 09:38:16 +0000
+++ pkg/dem/FlowEngine.hpp 2011-04-05 16:58:41 +0000
@@ -113,7 +113,6 @@
((double, Pressure_BACK_Boundary, 0,,"Pressure imposed on back boundary"))
((double, Pressure_LEFT_Boundary, 0,, "Pressure imposed on left boundary"))
((double, Pressure_RIGHT_Boundary, 0,, "Pressure imposed on right boundary"))
- ((int, id_sphere, 0,, "Average velocity will be computed for all cells incident to that sphere"))
((Vector3r, id_force, 0,, "Fluid force acting on sphere with id=flow.id_sphere"))
((bool, BOTTOM_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, TOP_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
