yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10162
[Branch ~yade-pkg/yade/git-trunk] Rev 3723: Save the relative velocities of particles.
------------------------------------------------------------
revno: 3723
committer: Donia <donia.marzougui@xxxxxxxxxxxxxxx>
timestamp: Mon 2013-10-21 15:11:33 +0200
message:
Save the relative velocities of particles.
modified:
lib/triangulation/FlowBoundingSphere.hpp
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.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/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2013-08-22 14:32:01 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2013-10-21 13:11:33 +0000
@@ -80,6 +80,10 @@
vector <Vector3r> normLubForce;
vector <Matrix3r> viscousBodyStress;
vector <Matrix3r> lubBodyStress;
+ vector <Vector3r> deltaNormVel;
+ vector <Vector3r> deltaShearVel;
+ vector <Vector3r> normalV;
+ vector <Real> surfaceDistance;
void Localize();
void Compute_Permeability();
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2013-09-23 17:39:59 +0000
+++ pkg/dem/FlowEngine.cpp 2013-10-21 13:11:33 +0000
@@ -87,7 +87,7 @@
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);
- if (viscousShear || shearLubrication){
+ if (shearLubrication || viscousShear){
force = force + solver->viscousShearForces[V_it->info().id()];
scene->forces.addTorque ( V_it->info().id(), solver->viscousShearTorques[V_it->info().id()]);
}
@@ -276,7 +276,7 @@
if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0)) flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
if ( WaveAction ) flow->ApplySinusoidalPressure ( flow->T[flow->currentTes].Triangulation(), sineMagnitude, sineAverage, 30 );
- if ( viscousShear || normalLubrication || shearLubrication) flow->computeEdgesSurfaces();
+ if (normalLubrication || shearLubrication) flow->computeEdgesSurfaces();
}
void FlowEngine::setPositionsBuffer(bool current)
@@ -582,7 +582,7 @@
template<class Solver>
void FlowEngine::ComputeViscousForces ( Solver& flow )
{
- if (viscousShear || normalLubrication || shearLubrication){
+ if (normalLubrication || shearLubrication || viscousShear){
if ( Debug ) cout << "Application of viscous forces" << endl;
if ( Debug ) cout << "Number of edges = " << flow.Edge_ids.size() << endl;
for ( unsigned int k=0; k<flow.viscousShearForces.size(); k++ ) flow.viscousShearForces[k]=Vector3r::Zero();
@@ -593,6 +593,7 @@
typedef typename Solver::Tesselation Tesselation;
const Tesselation& Tes = flow.T[flow.currentTes];
+ flow.deltaShearVel.clear(); flow.normalV.clear(); flow.deltaNormVel.clear(); flow.surfaceDistance.clear();
for ( int i=0; i< ( int ) flow.Edge_ids.size(); i++ ) {
const int& id1 = flow.Edge_ids[i].first;
@@ -600,6 +601,7 @@
int hasFictious= Tes.vertex ( id1 )->info().isFictious + Tes.vertex ( id2 )->info().isFictious;
if (hasFictious>0 or id1==id2) continue;
+ if (hasFictious==2) continue;
const shared_ptr<Body>& sph1 = Body::byId ( id1, scene );
const shared_ptr<Body>& sph2 = Body::byId ( id2, scene );
Sphere* s1=YADE_CAST<Sphere*> ( sph1->shape.get() );
@@ -648,26 +650,34 @@
}
}
deltaShearV = deltaV - ( normal.dot ( deltaV ) ) *normal;
+ flow.deltaShearVel.push_back(deltaShearV);
+ flow.normalV.push_back(normal);
+ flow.surfaceDistance.push_back(max(surfaceDist, 0.) + eps*meanRad);
+ if (affiche_vitesse) cout << "le vecteur normal entre " << id1 << " et " << id2 << "est " << normal <<" total " << ( int ) flow.Edge_ids.size()<< endl;
if (shearLubrication)
visc_f = flow.computeShearLubricationForce(deltaShearV,surfaceDist,i,eps,O1O2,meanRad);
else if (viscousShear)
visc_f = flow.computeViscousShearForce ( deltaShearV, i , Rh);
-
+
if (viscousShear || shearLubrication){
+
flow.viscousShearForces[id1]+=visc_f;
flow.viscousShearForces[id2]+=(-visc_f);
flow.viscousShearTorques[id1]+=O1C_vect.cross(visc_f);
flow.viscousShearTorques[id2]+=O2C_vect.cross(-visc_f);
+
/// Compute the viscous shear stress on each particle
if (viscousShearBodyStress){
flow.viscousBodyStress[id1] += visc_f * O1C_vect.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
flow.viscousBodyStress[id2] += (-visc_f) * O2C_vect.transpose()/ (4.0/3.0 *3.14* pow(r2,3));}
}
+
/// Compute the normal lubrication force applied on each particle
if (normalLubrication){
deltaNormV = normal.dot(deltaV);
+ flow.deltaNormVel.push_back(deltaNormV * normal);
lub_f = flow.computeNormalLubricationForce (deltaNormV, surfaceDist, i,eps,stiffness,scene->dt,meanRad)*normal;
flow.normLubForce[id1]+=lub_f;
flow.normLubForce[id2]+=(-lub_f);
@@ -677,7 +687,8 @@
flow.lubBodyStress[id1] += lub_f * O1C_vect.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
flow.lubBodyStress[id2] += (-lub_f) *O2C_vect.transpose() / (4.0/3.0 *3.14* pow(r2,3));}
}
-
+
+ if (affiche_force) cout<<"force tangentielle "<<visc_f<< " force normale "<< lub_f<<endl;
}
}
}
@@ -735,7 +746,7 @@
assert (Tes.vertexHandles[id] != NULL);
const Tesselation::Vertex_Info& v_info = Tes.vertexHandles[id]->info();
force =(pressureForce) ? Vector3r ( ( v_info.forces ) [0],v_info.forces[1],v_info.forces[2] ) : Vector3r(0,0,0);
- if (viscousShear){
+ if (shearLubrication || viscousShear){
force = force +solver->viscousShearForces[v_info.id()];
scene->forces.addTorque ( v_info.id(), solver->viscousShearTorques[v_info.id()]);
}
@@ -1089,7 +1100,7 @@
// if ( !first && (useSolver==0 || fluidBulkModulus>0)) flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
if ( WaveAction ) flow->ApplySinusoidalPressure ( Tes.Triangulation(), sineMagnitude, sineAverage, 30 );
- if ( viscousShear || normalLubrication || shearLubrication) flow->computeEdgesSurfaces();
+ if (normalLubrication || shearLubrication) flow->computeEdgesSurfaces();
if ( Debug ) cout << endl << "end buildTri------" << endl << endl;
}
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2013-08-30 13:07:40 +0000
+++ pkg/dem/FlowEngine.hpp 2013-10-21 13:11:33 +0000
@@ -95,6 +95,16 @@
return (flow->viscousBodyStress.size()>id_sph)?flow->viscousBodyStress[id_sph]:Matrix3r::Zero();}
TPL Matrix3r bodyNormalLubStress(unsigned int id_sph, Solver& flow) {
return (flow->lubBodyStress.size()>id_sph)?flow->lubBodyStress[id_sph]:Matrix3r::Zero();}
+ TPL Vector3r shearVelocity(unsigned int interaction, Solver& flow) {
+ return (flow->deltaShearVel[interaction]);}
+ TPL Vector3r normalVelocity(unsigned int interaction, Solver& flow) {
+ return (flow->deltaNormVel[interaction]);}
+ TPL Vector3r normalVect(unsigned int interaction, Solver& flow) {
+ return (flow->normalV[interaction]);}
+ TPL Real surfaceDistanceParticle(unsigned int interaction, Solver& flow) {
+ return (flow->surfaceDistance[interaction]);}
+ TPL Real edgeSize(Solver& flow) {
+ return (flow->Edge_ids.size());}
TPL python::list getConstrictions(bool all, Solver& flow) {
vector<Real> csd=flow->getConstrictions(); python::list pycsd;
for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k]!=0) || csd[k]>0) pycsd.append(csd[k]); return pycsd;}
@@ -147,6 +157,10 @@
Matrix3r _bodyShearLubStress(unsigned int id_sph) {return bodyShearLubStress(id_sph,solver);}
Matrix3r _bodyNormalLubStress(unsigned int id_sph) {return bodyNormalLubStress(id_sph,solver);}
Vector3r _fluidForce(unsigned int id_sph) {return fluidForce(id_sph,solver);}
+ Vector3r _shearVelocity(unsigned int interaction) {return shearVelocity(interaction,solver);}
+ Vector3r _normalVelocity(unsigned int interaction) {return normalVelocity(interaction,solver);}
+ Vector3r _normalVect(unsigned int interaction) {return normalVect(interaction,solver);}
+ Real _surfaceDistanceParticle(unsigned int interaction) {return surfaceDistanceParticle(interaction,solver);}
void _imposeFlux(Vector3r pos, Real flux) {return imposeFlux(pos,flux,*solver);}
unsigned int _imposePressure(Vector3r pos, Real p) {return imposePressure(pos,p,solver);}
void _setImposedPressure(unsigned int cond, Real p) {setImposedPressure(cond,p,solver);}
@@ -217,6 +231,9 @@
((int, ignoredBody,-1,,"Id of a sphere to exclude from the triangulation.)"))
((vector<int>, wallIds,vector<int>(6),,"body ids of the boundaries (default values are ok only if aabbWalls are appended before spheres, i.e. numbered 0,...,5)"))
((vector<bool>, boundaryUseMaxMin, vector<bool>(6,true),,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
+ ((bool, affiche_vitesse, false,,"show the relative velocities between particles"))
+ ((bool, affiche_force, false,,"show the lubrication force applied on particles"))
+ ((bool, create_file, false,,"create file of velocities"))
((bool, viscousShear, false,,"Compute viscous shear terms as developped by Donia Marzougui (FIXME: ref.)"))
((bool, shearLubrication, false,,"Compute shear lubrication force as developped by Brule (FIXME: ref.) "))
((double, eps, 0.00001,,"roughness defined as a fraction of particles size, giving the minimum distance between particles in the lubrication model."))
@@ -263,6 +280,11 @@
.def("normalLubForce",&FlowEngine::_normalLubForce,(python::arg("Id_sph")),"Return the normal lubrication force on sphere Id_sph.")
.def("bodyShearLubStress",&FlowEngine::_bodyShearLubStress,(python::arg("Id_sph")),"Return the shear lubrication stress on sphere Id_sph.")
.def("bodyNormalLubStress",&FlowEngine::_bodyNormalLubStress,(python::arg("Id_sph")),"Return the normal lubrication stress on sphere Id_sph.")
+ .def("shearVelocity",&FlowEngine::_shearVelocity,(python::arg("id_sph")),"Return the shear velocity of the interaction.")
+ .def("normalVelocity",&FlowEngine::_normalVelocity,(python::arg("id_sph")),"Return the normal velocity of the interaction.")
+ .def("normalVect",&FlowEngine::_normalVect,(python::arg("id_sph")),"Return the normal vector between particles.")
+ .def("surfaceDistanceParticle",&FlowEngine::_surfaceDistanceParticle,(python::arg("interaction")),"Return the distance between particles.")
+
.def("pressureProfile",&FlowEngine::pressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
.def("getPorePressure",&FlowEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
.def("averageSlicePressure",&FlowEngine::averageSlicePressure,(python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
@@ -343,6 +365,11 @@
Vector3r _normalLubForce(unsigned int id_sph) {return normalLubForce(id_sph,solver);}
Matrix3r _bodyShearLubStress(unsigned int id_sph) {return bodyShearLubStress(id_sph,solver);}
Matrix3r _bodyNormalLubStress(unsigned int id_sph) {return bodyNormalLubStress(id_sph,solver);}
+ Vector3r _shearVelocity(unsigned int id_sph) {return shearVelocity(id_sph,solver);}
+ Vector3r _normalVelocity(unsigned int id_sph) {return normalVelocity(id_sph,solver);}
+ Vector3r _normalVect(unsigned int id_sph) {return normalVect(id_sph,solver);}
+ Real _surfaceDistanceParticle(unsigned int interaction) {return surfaceDistanceParticle(interaction,solver);}
+ Real _edgeSize() {return edgeSize(solver);}
Vector3r _fluidForce(unsigned int id_sph) {return fluidForce(id_sph, solver);}
void _imposeFlux(Vector3r pos, Real flux) {return imposeFlux(pos,flux,*solver);}
@@ -395,6 +422,11 @@
.def("normalLubForce",&PeriodicFlowEngine::_normalLubForce,(python::arg("Id_sph")),"Return the normal lubrication force on sphere Id_sph.")
.def("bodyShearLubStress",&PeriodicFlowEngine::_bodyShearLubStress,(python::arg("Id_sph")),"Return the shear lubrication stress on sphere Id_sph.")
.def("bodyNormalLubStress",&PeriodicFlowEngine::_bodyNormalLubStress,(python::arg("Id_sph")),"Return the normal lubrication stress on sphere Id_sph.")
+ .def("shearVelocity",&PeriodicFlowEngine::_shearVelocity,(python::arg("id_sph")),"Return the shear velocity of the interaction.")
+ .def("normalVelocity",&PeriodicFlowEngine::_normalVelocity,(python::arg("id_sph")),"Return the normal velocity of the interaction.")
+ .def("normalVect",&PeriodicFlowEngine::_normalVect,(python::arg("id_sph")),"Return the normal vector between particles.")
+ .def("surfaceDistanceParticle",&PeriodicFlowEngine::_surfaceDistanceParticle,(python::arg("interaction")),"Return the distance between particles.")
+ .def("edgeSize",&PeriodicFlowEngine::_edgeSize,"Return the number of interactions.")
// .def("imposeFlux",&FlowEngine::_imposeFlux,(python::arg("pos"),python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
.def("saveVtk",&PeriodicFlowEngine::saveVtk,"Save pressure field in vtk format.")
