yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11381
[Branch ~yade-pkg/yade/git-trunk] Rev 3377: Merge github.com:yade/trunk into chaoUnsat
Merge authors:
Anton Gladky (gladky-anton)
Bruno Chareyre (bruno-chareyre)
Christian Jakob (jakob-ifgt)
Jan Stránský (honzik)
Kneib François (francois-kneib)
------------------------------------------------------------
revno: 3377 [merge]
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Wed 2013-06-19 14:47:52 +0200
message:
Merge github.com:yade/trunk into chaoUnsat
added:
examples/clumps/apply-buoyancy-clumps.py
examples/rotationalResistance.py
pkg/dem/BubbleMat.cpp
pkg/dem/BubbleMat.hpp
pkg/dem/ViscoelasticCapillarPM.cpp
modified:
doc/references.bib
doc/sphinx/references.bib
doc/yade-articles.bib
examples/CapillaryPhys-example.py
lib/triangulation/FlowBoundingSphere.ipp
pkg/common/Grid.cpp
pkg/common/Grid.hpp
pkg/common/InsertionSortCollider.cpp
pkg/common/OpenGLRenderer.cpp
pkg/common/OpenGLRenderer.hpp
pkg/dem/CapillaryPhys.hpp
pkg/dem/CapillaryStressRecorder.cpp
pkg/dem/CapillaryTriaxialTest.cpp
pkg/dem/CapillaryTriaxialTest.hpp
pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp
pkg/dem/SampleCapillaryPressureEngine.cpp
pkg/dem/Shop.cpp
pkg/dem/Shop.hpp
pkg/dem/ViscoelasticPM.cpp
pkg/dem/ViscoelasticPM.hpp
py/_utils.cpp
--
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 'doc/references.bib'
--- doc/references.bib 2013-05-28 06:22:52 +0000
+++ doc/references.bib 2013-06-12 14:35:49 +0000
@@ -571,3 +571,19 @@
URL = {http://www.tandfonline.com/doi/abs/10.1080/00018730500167855},
eprint = {http://www.tandfonline.com/doi/pdf/10.1080/00018730500167855}
}
+
+@article{Rabinov2005,
+ author = "RABINOVICH Yakov I. and ESAYANUR Madhavan S. and MOUDGIL Brij M.",
+ institution = "Particle Engineering Research Center, University of Florida, USA; Department of Materials Science and Engineering, University of Florida, USA",
+ title = "Capillary forces between two spheres with a fixed volume liquid bridge : Theory and experiment",
+ journal = "Langmuir",
+ year = "2005",
+ volume = "21",
+ number = "24",
+ pages = "10992--10997",
+ editor = "American Chemical Society",
+ issn = "0743-7463",
+ note = "eng",
+ keywords = "Atomic force microscopy; Equation; Plane surface; Rupture; Thermodynamics; Methodology; Gas pressure; Calculation; Solid; Interaction energy; Geometry; Plate; Separation; Prediction; Powder; Vapor; Capillary condensation; Theory; Liquid; Force",
+ url = "http://www.refdoc.fr/Detailnotice?idarticle=7435486";
+}
=== modified file 'doc/sphinx/references.bib'
--- doc/sphinx/references.bib 2010-06-12 20:39:00 +0000
+++ doc/sphinx/references.bib 2013-06-06 09:27:41 +0000
@@ -269,7 +269,6 @@
Month = {October},
}
-
@inproceedings{Price2007,
Author = { Mathew Price and Vasile Murariu and Garry Morrison},
title= {Sphere clump generation and trajectory comparison for real particles},
@@ -279,19 +278,39 @@
}
@inproceedings{Kuhl2001,
- title = {Microplane modelling and particle modelling of cohesive-frictional materials},
- author = {E. Kuhl and G. A. D'Addetta and M. Leukart and E. Ramm},
- booktitle = {Continuous and Discontinuous Modelling of Cohesive-Frictional Materials},
- publisher = {Springer Berlin / Heidelberg},
- issn = {1616-6361},
- isbn = {978-3-540-41525-1},
- url = {http://www.springerlink.com/content/e50544266r506615},
- abstract = {This paper aims at comparing the microplane model as a particular representative of the class of continuous material models with a discrete particle model which is of discontinuous nature. Thereby, the constitutive equations of both approaches will be based onVoigtshypothesis defining the strain state on the individual microplanes as well as the relative particle displacements. Through an appropriate constitutive assumption, the microplane stresses and the contact forces can be determined. In both cases, the equivalence of microscopic and macroscopic virtual work yields the overall stress strain relation. An elastic and an elasto-plastic material characterization of the microplane model and the particle model are derived and similarities of both approaches are illustrated.},
- volume = {568},
- series = {Lecture Notes in Physics},
- year = {2001},
- pages = {31--46},
- doi = {10.1007/3-540-44424-6_3},
+ title = {Microplane modelling and particle modelling of cohesive-frictional materials},
+ author = {E. Kuhl and G. A. D'Addetta and M. Leukart and E. Ramm},
+ booktitle = {Continuous and Discontinuous Modelling of Cohesive-Frictional Materials},
+ publisher = {Springer Berlin / Heidelberg},
+ issn = {1616-6361},
+ isbn = {978-3-540-41525-1},
+ url = {http://www.springerlink.com/content/e50544266r506615},
+ abstract = {This paper aims at comparing the microplane model as a particular representative of the class of continuous material models with a discrete particle model which is of discontinuous nature. Thereby, the constitutive equations of both approaches will be based onVoigtshypothesis defining the strain state on the individual microplanes as well as the relative particle displacements. Through an appropriate constitutive assumption, the microplane stresses and the contact forces can be determined. In both cases, the equivalence of microscopic and macroscopic virtual work yields the overall stress strain relation. An elastic and an elasto-plastic material characterization of the microplane model and the particle model are derived and similarities of both approaches are illustrated.},
+ volume = {568},
+ series = {Lecture Notes in Physics},
+ year = {2001},
+ pages = {31--46},
+ doi = {10.1007/3-540-44424-6_3},
}
+@article{Zhou1999536,
+ title = "Rolling friction in the dynamic simulation of sandpile formation ",
+ journal = "Physica A: Statistical Mechanics and its Applications ",
+ volume = "269",
+ number = "2–4",
+ pages = "536 - 553",
+ year = "1999",
+ note = "",
+ issn = "0378-4371",
+ doi = "10.1016/S0378-4371(99)00183-1",
+ url = "http://www.sciencedirect.com/science/article/pii/S0378437199001831";,
+ author = "Y.C. Zhou and B.D. Wright and R.Y. Yang and B.H. Xu and A.B. Yu",
+ keywords = "Static sandpiles",
+ keywords = "Granular compaction",
+ keywords = "Dynamics and kinematics of a particle and a system of particles",
+ keywords = "Granular flow",
+ keywords = "Mixing",
+ keywords = "Segregation and stratification",
+ keywords = "Granular systems "
+}
=== modified file 'doc/yade-articles.bib'
--- doc/yade-articles.bib 2013-05-03 16:45:52 +0000
+++ doc/yade-articles.bib 2013-06-12 14:30:56 +0000
@@ -485,3 +485,21 @@
author = {François Nicot and Nejib Hadda and Mohamed Guessasma and Jerome Fortin and Olivier Millet}
}
+@article{Gusenbauer2012,
+ title={Self-organizing magnetic beads for biomedical applications},
+ author={Gusenbauer, Markus and Kovacs, Alexander and Reichel, Franz and Exl, Lukas and Bance, Simon and {\"O}zelt, Harald and Schrefl, Thomas},
+ journal={Journal of Magnetism and Magnetic Materials},
+ volume={324},
+ number={6},
+ pages={977--982},
+ year={2012},
+ publisher={Elsevier}
+}
+
+@article{Bourrier2013,
+ title = {Discrete modeling of granular soils reinforcement by plant roots},
+ journal = {Ecological Engineering},
+ year = {2013},
+ doi = {10.1016/j.ecoleng.2013.05.002},
+ url = {http://dx.doi.org/10.1016/j.ecoleng.2013.05.002}
+}
=== modified file 'examples/CapillaryPhys-example.py'
--- examples/CapillaryPhys-example.py 2013-03-28 10:53:06 +0000
+++ examples/CapillaryPhys-example.py 2013-05-30 20:17:45 +0000
@@ -45,7 +45,7 @@
[Ip2_FrictMat_FrictMat_MindlinCapillaryPhys(label='ContactModel')],#for hertz model only
[Law2_ScGeom_MindlinPhys_Mindlin()]#for hertz model only
),
- Law2_ScGeom_CapillaryPhys_Capillarity(CapillaryPressure=10000),#for hertz model only
+ Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=10000),#for hertz model only
NewtonIntegrator(damping=local_damping,gravity=(0,0,-9.81)),
]
ContactModel.betan=viscous_normal
@@ -60,7 +60,7 @@
[Ip2_FrictMat_FrictMat_CapillaryPhys()], #for linear model only
[Law2_ScGeom_FrictPhys_CundallStrack()], #for linear model only
),
- Law2_ScGeom_CapillaryPhys_Capillarity(CapillaryPressure=10000),#for linear model only
+ Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=10000),#for linear model only
NewtonIntegrator(damping=local_damping,gravity=(0,0,-9.81)),
]
=== added file 'examples/clumps/apply-buoyancy-clumps.py'
--- examples/clumps/apply-buoyancy-clumps.py 1970-01-01 00:00:00 +0000
+++ examples/clumps/apply-buoyancy-clumps.py 2013-06-15 08:19:25 +0000
@@ -0,0 +1,145 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+''' The script shows how to include the effect of buoyancy in a particle assembly
+ under condition of an increasing water-level. Water-level at start of the sim-
+ ulation is at the lower boundary of the model. During the calculation particles
+ gets partly surrounded by water and buoyancy (=density*volumeOfDisplacedWater)
+ is increasing until particle is completely surrounded. When particle is sur-
+ rounded volumeOfDisplacedWater is equal to the volume of the particle.
+
+ For calculation of buoyancy of a clump the theoretical radius
+ R = (3*mass/(4*pi*density))^1/3
+ of a sphere with clump mass and clump volume
+ V = (4*pi*R^3)/3 = mass/density
+ is used. This approximation can be used for well rounded clumps.
+
+ Buoyancy is included via reduced mass (=massAtStart - dm), where dm is the
+ buoyancy. Reduced mass must be corrected via density scaling for calculation
+ of correct particle accelerations.'''
+
+#define material properties:
+shear_modulus = 3.2e10
+poisson_ratio = 0.15
+young_modulus = 2*shear_modulus*(1+poisson_ratio)
+angle = 0.5 #friction angle in radians
+rho_p = 2650 #density of particles
+
+v_iw = 1 #velocity of increasing water-level
+
+#model boundaries:
+boundaryMin = Vector3(-1.5,-1.5,0)
+boundaryMax = Vector3(1.5,1.5,2)
+
+#material:
+id_Mat=O.materials.append(FrictMat(young=young_modulus,poisson=poisson_ratio,density=rho_p,frictionAngle=angle))
+Mat=O.materials[id_Mat]
+
+#create assembly of spheres:
+sp=pack.SpherePack()
+sp.makeCloud(minCorner=boundaryMin,maxCorner=boundaryMax,rMean=.2,rRelFuzz=.5,num=100,periodic=False)
+O.bodies.append([sphere(c,r,material=Mat) for c,r in sp])
+
+#create template for clumps and replace 10% of spheres by clumps:
+templateList = [clumpTemplate(relRadii=[1,.5],relPositions=[[0,0,0],[.7,0,0]])]
+O.bodies.replaceByClumps(templateList,[.1])
+
+#create boundary:
+O.bodies.append(facetBox((0,0,1), (boundaryMax-boundaryMin)/2, fixed=True, material=Mat))#boundaryMax-boundaryMin
+
+#define engines:
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
+ [Ip2_FrictMat_FrictMat_MindlinPhys(betan=0.0,betas=0.0,label='ContactModel')],
+ [Law2_ScGeom_MindlinPhys_Mindlin(neverErase=False,label='Mindlin')]
+ ),
+ GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8, defaultDt=0.9*PWaveTimeStep(),label='ts'),
+ NewtonIntegrator(gravity=(0,0,-10),damping=0.7,label='integrator')
+]
+ts.dead=True
+
+#definition to apply buoyancy:
+def apply_buo(water_height,saturatedList,startMass):
+ for b in O.bodies:
+ if b not in saturatedList:
+ h_low = 1e9
+ h_high = -1e9
+ if b.isStandalone and isinstance(b.shape,Sphere):
+ rad = b.shape.radius
+ pos = b.state.pos
+ h_low = pos[2] - rad
+ h_high = pos[2] + rad
+ elif b.isClump: #determine rad, h_high and h_low for clumps:
+ keys = O.bodies[b.id].shape.members.keys()
+ for ii in range(0,len(keys)):
+ pos = O.bodies[keys[ii]].state.pos
+ h_low = min(h_low,pos[2]-O.bodies[keys[ii]].shape.radius)
+ h_high = max(h_high,pos[2]+O.bodies[keys[ii]].shape.radius)
+ rad = ( 3*startMass[b.id]/(4*pi*O.bodies[keys[ii]].mat.density) )**(1./3.) #get radius from startMass
+ else:
+ continue
+ if water_height > h_low:
+ dh = water_height - h_low
+ dh = min(dh,2*rad) #to get sure, that dh is not bigger than 2*radius
+ dm = (pi/3)*dh*dh*(3*rad - dh)*1000 #buoyancy acts via reduced mass dm=rho*V_cap
+ b.state.mass = startMass[b.id] - dm #needs a list of masses of all particles at start (see 5-final.py)
+ #reduced mass must be corrected via density scaling for accelerations a = F/m_reduced, a_correct = a*densityScaling, see line 179 in pkg/dem/NewtonIntegrator.cpp:
+ b.state.densityScaling = (startMass[b.id] - dm)/startMass[b.id]
+ if water_height > h_high:
+ saturatedList.append(b)
+ return saturatedList
+
+#STEP1: reduce overlaps from replaceByClumps() method:
+O.dt=1e-6 #small time step for preparation steps via calm()
+print '\nSTEP1 in progress. Please wait a minute ...\n'
+O.engines=O.engines+[PyRunner(iterPeriod=10000,command='calm()',label='calmRunner')]
+O.run(1000000,True)
+calmRunner.dead=True
+
+#STEP2: let particles settle down:
+O.dt=1e-5
+print '\nSTEP2 in progress. Please wait a minute ...\n'
+O.run(50000,True)
+
+#get maximum body id:
+idmax=0
+for b in O.bodies:
+ idmax=max(idmax,b.id)
+integrator.densityScaling=True #activate density scaling
+
+#get a list of masses at start:
+startMass = [0 for ii in range(0,idmax+1)]
+for b in O.bodies:
+ if (b.isStandalone and isinstance(b.shape,Sphere)) or b.isClump:
+ startMass[b.id] = b.state.mass
+
+#start PyRunner engine to apply buoyancy:
+saturatedList = []
+t0 = O.time
+O.engines=O.engines+[PyRunner(iterPeriod=100,command='saturatedList=apply_buo(((O.time-t0) * v_iw),saturatedList,startMass)',label='buoLabel')]
+
+#create 2 facets, that show water height:
+idsWaterFacets = []
+idsWaterFacets.append(O.bodies.append(facet( \
+ [ boundaryMin, [boundaryMax[0],boundaryMin[1],boundaryMin[2]], [boundaryMax[0],boundaryMax[1],boundaryMin[2]] ], \
+ fixed=True,material=FrictMat(young=0),color=(0,0,1),wire=False)))#no interactions will appear
+idsWaterFacets.append(O.bodies.append(facet( \
+ [ [boundaryMax[0],boundaryMax[1],boundaryMin[2]], [boundaryMin[0],boundaryMax[1],boundaryMin[2]], boundaryMin ], \
+ fixed=True,material=FrictMat(young=0),color=(0,0,1),wire=False)))#no interactions will appear
+
+#set velocity of incr. water level to the facets:
+for ids in idsWaterFacets:
+ O.bodies[ids].state.vel = [0,0,v_iw]
+
+#STEP3: simulate buoyancy#
+O.dt=3e-5
+from yade import qt
+qt.Controller()
+v=qt.View()
+v.eyePosition=(-7,0,2); v.upVector=(0,0,1); v.viewDir=(1,0,-.1); v.axes=True; v.sceneRadius=1.9
+print '\nSTEP3 started ...\n'
+O.run(60000)
+
=== added file 'examples/rotationalResistance.py'
--- examples/rotationalResistance.py 1970-01-01 00:00:00 +0000
+++ examples/rotationalResistance.py 2013-06-06 09:27:41 +0000
@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+# encoding: utf-8
+from yade import utils, plot
+o = Omega()
+fr = 0.5;rho=2000
+tc = 0.001; en = 0.7; et = 0.7; o.dt = 0.0002*tc
+
+r = 0.002
+
+param = getViscoelasticFromSpheresInteraction(tc,en,et)
+
+mat1 = O.materials.append(ViscElMat(frictionAngle=fr,mR = 0.05, mRtype = 1, density=rho,**param))
+mat2 = O.materials.append(ViscElMat(frictionAngle=fr,mR = 0.05, mRtype = 2, density=rho,**param))
+
+oriBody = Quaternion(Vector3(1,0,0),(pi/28))
+
+id1 = O.bodies.append(sphere(center=[0,0,2*r],radius=r,material=mat1))
+id2 = O.bodies.append(geom.facetBox(center=(0,-16.0*r,-2*r),orientation=oriBody,extents=(r,17.0*r,0), material=mat1,color=(0,0,1)))
+
+id3 = O.bodies.append(sphere(center=[10*r,0,2*r],radius=r,material=mat2))
+id4 = O.bodies.append(geom.facetBox(center=(10*r,-16.0*r,-2*r),orientation=oriBody,extents=(r,17.0*r,0), material=mat2,color=(0,0,1)))
+
+o.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
+ [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
+ [Law2_ScGeom_ViscElPhys_Basic()],
+ ),
+ NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),
+ PyRunner(command='addPlotData()',iterPeriod=10000, dead = False, label='graph'),
+]
+
+def addPlotData():
+ f1 = [0,0,0]
+ s1 = O.bodies[id1].state.pos[1]
+ s2 = O.bodies[id3].state.pos[1]
+ plot.addData(sc=O.time, fc1=s1, fc2=s2)
+
+
+
+plot.plots={'sc':('fc1','fc2')}; plot.plot()
+
+from yade import qt
+qt.View()
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2013-05-03 16:45:52 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2013-06-07 07:47:59 +0000
@@ -897,8 +897,13 @@
{
RTriangulation& Tri = T[currentTes].Triangulation();
vector<double> constrictions;
- for (Finite_facets_iterator f_it=Tri.finite_facets_begin(); f_it != Tri.finite_facets_end();f_it++)
+ for (Finite_facets_iterator f_it=Tri.finite_facets_begin(); f_it != Tri.finite_facets_end();f_it++){
+ //in the periodic case, we skip facets with lowest id out of the base period
+ if ( ((f_it->first->info().index < f_it->first->neighbor(f_it->second)->info().index) && f_it->first->info().isGhost)
+ || ((f_it->first->info().index > f_it->first->neighbor(f_it->second)->info().index) && f_it->first->neighbor(f_it->second)->info().isGhost)
+ || f_it->first->info().index == 0 || f_it->first->neighbor(f_it->second)->info().index == 0) continue;
constrictions.push_back(Compute_EffectiveRadius(f_it->first, f_it->second));
+ }
return constrictions;
}
@@ -908,6 +913,10 @@
RTriangulation& Tri = T[currentTes].Triangulation();
vector<Constriction> constrictions;
for (Finite_facets_iterator f_it=Tri.finite_facets_begin(); f_it != Tri.finite_facets_end();f_it++){
+ //in the periodic case, we skip facets with lowest id out of the base period
+ if ( ((f_it->first->info().index < f_it->first->neighbor(f_it->second)->info().index) && f_it->first->info().isGhost)
+ || ((f_it->first->info().index > f_it->first->neighbor(f_it->second)->info().index) && f_it->first->neighbor(f_it->second)->info().isGhost)
+ || f_it->first->info().index == 0 || f_it->first->neighbor(f_it->second)->info().index == 0) continue;
vector<double> rn;
const Vecteur& normal = f_it->first->info().facetSurfaces[f_it->second];
if (!normal[0] && !normal[1] && !normal[2]) continue;
=== modified file 'pkg/common/Grid.cpp'
--- pkg/common/Grid.cpp 2013-04-30 09:52:44 +0000
+++ pkg/common/Grid.cpp 2013-06-03 15:28:41 +0000
@@ -24,6 +24,9 @@
ScGridCoGeom::~ScGridCoGeom(){}
YADE_PLUGIN((ScGridCoGeom));
+GridCoGridCoGeom::~GridCoGridCoGeom(){}
+YADE_PLUGIN((GridCoGridCoGeom));
+
void GridNode::addConnection(shared_ptr<Body> GC){
ConnList.push_back(GC);
}
@@ -69,6 +72,96 @@
}
YADE_PLUGIN((Ig2_GridNode_GridNode_GridNodeGeom6D));
+//! \\//
+bool Ig2_GridConnection_GridConnection_GridCoGridCoGeom::go( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c)
+{
+ /*FIXME : /!\ Note that this geometry doesn't take care of any unwished duplicated contact or shear force following. /!\*/
+ GridConnection* conn1 = YADE_CAST<GridConnection*>(cm1.get());
+ GridConnection* conn2 = YADE_CAST<GridConnection*>(cm2.get());
+ State* stNode11 = conn1->node1->state.get();
+ State* stNode12 = conn1->node2->state.get();
+ State* stNode21 = conn2->node1->state.get();
+ State* stNode22 = conn2->node2->state.get();
+
+ if(conn1->node1==conn2->node1 || conn1->node1==conn2->node2 || conn1->node2==conn2->node1 || conn1->node2==conn2->node2){
+ //Two connections share at least one node, so they are contiguous => they must not interact.
+ return false;
+ }
+ //There could be a contact between to connections. Check this now.
+ bool isNew = !c->geom;
+ Real k,m;
+ Vector3r A=stNode11->pos, a=stNode12->pos-A; //"A" is an extremity of conn1, "a" is the connection's segment.
+ Vector3r B=stNode21->pos, b=stNode22->pos-B; //"B" is an extremity of conn2, "b" is the connection's segment.
+ B+=shift2;//periodicity.
+ Vector3r N=a.cross(b); //"N" is orthogonal to "a" and "b". It means that "N" describes the common plan between a and b.
+ if(N.norm()>1e-14){ //If "a" and "b" are colinear, "N==0" and this is a special case.
+ Real dist=N.dot(B-A)/(N.norm()); //here "dist" is oriented, so it's sign depends on the orientation of "N" against "AB".
+ Vector3r pB=B-dist*(N/(N.norm())); //"pB" is the projection of the point "B" in the plane defined by his normal vector "N".
+ //Now we have pB, so we will compute the intersection of two segments into a plane.
+ int b0, b1; //2 base vectors used to compute the segment intersection. For more accuracy and to avoid det==0, don't choose the axis where N is max.
+ if(abs(N[0])<abs(N[1]) || abs(N[0])<abs(N[2])){b0=0 ; b1=abs(N[1])<abs(N[2])?1:2;}
+ else { b0=1;b1=2;}
+ Real det=a[b0]*b[b1]-a[b1]*b[b0];
+ if (abs(det)>1e-14){
+ //Now compute k and m, who are the parameters (relative position on the connections) of the intersection on conn1 ("A" and "a") and conn2 ("B" and "b") respectively.
+ k = (b[b1]*(pB[b0]-A[b0])+b[b0]*(A[b1]-pB[b1]))/det;
+ m = (a[b0]*(-pB[b1]+A[b1])+a[b1]*(pB[b0]-A[b0]))/det;
+ //This is a little bit tricky : if we haven't 0<k,m<1, it means that the intersection is not inside both segments,
+ //but the contact can occurs anyway between a connection's extremity and a connection's edge or between two connection's extremity.
+ //So the three next lines : don't modify k and m if (0<k,m<1), but modify them otherwise to compute later the right normal and penetrationDepth of the contact.
+ k = max(min( k,1.0),0.0);
+ m = max(min( (A+a*k-B).dot(b)/(pow(b.norm(),2.0)) ,1.0),0.0);
+ k = max(min( (B+b*m-A).dot(a)/(pow(a.norm(),2.0)) ,1.0),0.0);
+ }
+ else {//should never happen
+ k=0;m=0;
+ cout<<"Error in Ig2_GridConnection_GridConnection_GridCoGridCoGeom : det=="<<det<<endl;
+ cout<<"Details : N="<<N<<" b0="<<b0<<" b1="<<b1<<" a="<<a<<" b="<<b<<endl;
+ }
+ }
+ else{ //this is a special case for perfectly colinear vectors ("a" and "b")
+ Real PA=(A-B).dot(b)/(b.norm()*b.norm()); PA=min(1.0,max(0.0,PA));
+ Real Pa=(A+a-B).dot(b)/(b.norm()*b.norm()); Pa=min(1.0,max(0.0,Pa));
+ Real PB=(B-A).dot(a)/(a.norm()*a.norm()); PB=min(1.0,max(0.0,PB));
+ Real Pb=(B+b-A).dot(a)/(a.norm()*a.norm()); Pb=min(1.0,max(0.0,Pb));
+ k=(PB+Pb)/2. ; m=(PA+Pa)/2.;
+ }
+
+ //Compute the geometry if "penetrationDepth" is positive.
+ double penetrationDepth = conn1->radius + conn2->radius - (A+k*a - (B+m*b)).norm();
+ shared_ptr<GridCoGridCoGeom> scm;
+ if(isNew){
+ if(penetrationDepth<0)return false;
+ scm=shared_ptr<GridCoGridCoGeom>(new GridCoGridCoGeom());
+ c->geom=scm;
+ }
+ else scm=YADE_PTR_CAST<GridCoGridCoGeom>(c->geom);
+ //k and m are used to compute almost everything...
+ //Fictious states (spheres) are generated at k or m of each connection, they will handle the contact.
+ scm->relPos1=k ; scm->relPos2=m;
+ scm->fictiousState1.pos=A + k*a ; scm->fictiousState2.pos=B + m*b;
+ scm->radius1 = conn1->radius ; scm->radius2 = conn2->radius;
+ scm->fictiousState1.vel = (1-k)*stNode11->vel + k*stNode12->vel;
+ scm->fictiousState2.vel = (1-m)*stNode21->vel + m*stNode22->vel;
+ Vector3r direction = a/(a.norm());
+ scm->fictiousState1.angVel = ((1-k)*stNode11->angVel + k*stNode12->angVel).dot(direction)*direction //twist part : interpolated
+ + a.cross(stNode12->vel - stNode11->vel);// non-twist part : defined from nodes velocities
+ direction = b/(b.norm());
+ scm->fictiousState2.angVel = ((1-m)*stNode21->angVel + m*stNode22->angVel).dot(direction)*direction //twist part : interpolated
+ + b.cross(stNode22->vel - stNode21->vel);// non-twist part : defined from nodes velocities
+ Vector3r normal= scm->fictiousState2.pos - scm->fictiousState1.pos;
+ normal/=normal.norm();
+ scm->contactPoint = scm->fictiousState1.pos + (scm->radius1-0.5*penetrationDepth)*normal;
+ scm->penetrationDepth=penetrationDepth;
+ scm->precompute(scm->fictiousState1,scm->fictiousState2,scene,c,normal,isNew,shift2,true);
+ return true;
+}
+
+bool Ig2_GridConnection_GridConnection_GridCoGridCoGeom::goReverse( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c)
+{
+ return go(cm1,cm2,state2,state1,-shift2,force,c);
+}
+YADE_PLUGIN((Ig2_GridConnection_GridConnection_GridCoGridCoGeom));
//! O/
bool Ig2_Sphere_GridConnection_ScGridCoGeom::go( const shared_ptr<Shape>& cm1,
@@ -403,6 +496,63 @@
}
YADE_PLUGIN((Law2_ScGridCoGeom_CohFrictPhys_CundallStrack));
+void Law2_GridCoGridCoGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+ int id1 = contact->getId1(), id2 = contact->getId2();
+ id_t id11 = (static_cast<GridConnection*>((&Body::byId(id1)->shape)->get()))->node1->getId();
+ id_t id12 = (static_cast<GridConnection*>((&Body::byId(id1)->shape)->get()))->node2->getId();
+ id_t id21 = (static_cast<GridConnection*>((&Body::byId(id2)->shape)->get()))->node1->getId();
+ id_t id22 = (static_cast<GridConnection*>((&Body::byId(id2)->shape)->get()))->node2->getId();
+ GridCoGridCoGeom* geom= static_cast<GridCoGridCoGeom*>(ig.get());
+ FrictPhys* phys = static_cast<FrictPhys*>(ip.get());
+ if(geom->penetrationDepth <0){
+ if (neverErase) {
+ phys->shearForce = Vector3r::Zero();
+ phys->normalForce = Vector3r::Zero();}
+ else scene->interactions->requestErase(contact);
+ return;}
+ Real& un=geom->penetrationDepth;
+ phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
+
+ Vector3r& shearForce = geom->rotate(phys->shearForce);
+ const Vector3r& shearDisp = geom->shearIncrement();
+ shearForce -= phys->ks*shearDisp;
+ Real maxFs = phys->normalForce.squaredNorm()*std::pow(phys->tangensOfFrictionAngle,2);
+
+ if (!scene->trackEnergy && !traceEnergy){//Update force but don't compute energy terms (see below))
+ // PFC3d SlipModel, is using friction angle. CoulombCriterion
+ if( shearForce.squaredNorm() > maxFs ){
+ Real ratio = sqrt(maxFs) / shearForce.norm();
+ shearForce *= ratio;}
+ } else {
+ //almost the same with additional Vector3r instatinated for energy tracing,
+ //duplicated block to make sure there is no cost for the instanciation of the vector when traceEnergy==false
+ if(shearForce.squaredNorm() > maxFs){
+ Real ratio = sqrt(maxFs) / shearForce.norm();
+ Vector3r trialForce=shearForce;//store prev force for definition of plastic slip
+ //define the plastic work input and increment the total plastic energy dissipated
+ shearForce *= ratio;
+ Real dissip=((1/phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/;
+ if (traceEnergy) plasticDissipation += dissip;
+ else if(dissip>0) scene->energy->add(dissip,"plastDissip",plastDissipIx,/*reset*/false);
+ }
+ // compute elastic energy as well
+ scene->energy->add(0.5*(phys->normalForce.squaredNorm()/phys->kn+phys->shearForce.squaredNorm()/phys->ks),"elastPotential",elastPotentialIx,/*reset at every timestep*/true);
+ }
+ Vector3r force = -phys->normalForce-shearForce;
+ Vector3r torque1 = (geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force);
+ Vector3r torque2 = (geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force);
+
+ scene->forces.addForce(id11,(1-geom->relPos1)*force);
+ scene->forces.addForce(id12,geom->relPos1*force);
+ scene->forces.addForce(id21,-(1-geom->relPos2)*force);
+ scene->forces.addForce(id22,-geom->relPos2*force);
+
+ scene->forces.addTorque(id11,(1-geom->relPos1)*torque1);
+ scene->forces.addTorque(id12,geom->relPos1*torque1);
+ scene->forces.addTorque(id21,(1-geom->relPos2)*torque2);
+ scene->forces.addTorque(id22,geom->relPos2*torque2);
+}
+YADE_PLUGIN((Law2_GridCoGridCoGeom_FrictPhys_CundallStrack));
//!################## Bounds #####################
void Bo1_GridConnection_Aabb::go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r& se3, const Body* b){
=== modified file 'pkg/common/Grid.hpp'
--- pkg/common/Grid.hpp 2013-05-14 21:24:11 +0000
+++ pkg/common/Grid.hpp 2013-06-03 15:28:41 +0000
@@ -28,6 +28,7 @@
#include <yade/core/Body.hpp>
#include <yade/pkg/common/Dispatching.hpp>
#include <yade/pkg/dem/Ig2_Sphere_Sphere_ScGeom.hpp>
+#include <yade/pkg/dem/ElasticContactLaw.hpp>
#ifdef YADE_OPENGL
#include<yade/pkg/common/GLDrawFunctors.hpp>
#endif
@@ -103,6 +104,21 @@
};
REGISTER_SERIALIZABLE(ScGridCoGeom);
+//! -|-
+class GridCoGridCoGeom: public ScGeom {
+ public:
+ /// Emulate a sphere whose position is the projection of sphere's center on cylinder sphere, and with motion linearly interpolated between nodes
+ State fictiousState1,fictiousState2;
+ virtual ~GridCoGridCoGeom ();
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(GridCoGridCoGeom,ScGeom,"Geometry of a :yref:`GridConnection`-:yref:`GridConnection` contact.",
+ ((Real,relPos1,0,,"position of the contact on the first connection (0: node-, 1:node+) |yupdate|"))
+ ((Real,relPos2,0,,"position of the contact on the first connection (0: node-, 1:node+) |yupdate|")),
+ createIndex(); /*ctor*/
+ );
+ REGISTER_CLASS_INDEX(ScGridCoGeom,ScGeom);
+};
+REGISTER_SERIALIZABLE(GridCoGridCoGeom);
+
//!################## IGeom Functors #####################
//! O-O
@@ -120,6 +136,19 @@
};
REGISTER_SERIALIZABLE(Ig2_GridNode_GridNode_GridNodeGeom6D);
+//! -/-
+class Ig2_GridConnection_GridConnection_GridCoGridCoGeom: public IGeomFunctor{
+ public:
+ virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ virtual bool goReverse( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ YADE_CLASS_BASE_DOC_ATTRS(Ig2_GridConnection_GridConnection_GridCoGridCoGeom,IGeomFunctor,"Create/update a :yref:`GridCoGridCoGeom` instance representing the geometry of a contact point between two :yref:`GridConnection` , including relative rotations.",
+ );
+ FUNCTOR2D(GridConnection,GridConnection);
+ DEFINE_FUNCTOR_ORDER_2D(GridConnection,GridConnection);
+};
+REGISTER_SERIALIZABLE(Ig2_GridConnection_GridConnection_GridCoGridCoGeom);
+
+
//! O/
class Ig2_Sphere_GridConnection_ScGridCoGeom: public IGeomFunctor{
public:
@@ -162,6 +191,18 @@
};
REGISTER_SERIALIZABLE(Law2_ScGridCoGeom_CohFrictPhys_CundallStrack);
+//! -/-
+class Law2_GridCoGridCoGeom_FrictPhys_CundallStrack: public Law2_ScGeom_FrictPhys_CundallStrack{
+ public:
+ virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ YADE_CLASS_BASE_DOC_ATTRS(Law2_GridCoGridCoGeom_FrictPhys_CundallStrack,Law2_ScGeom_FrictPhys_CundallStrack,"Frictional elastic contact law between two :yref:`gridConnection` . See :yref:`Law2_ScGeom_FrictPhys_CundallStrack` for more details.",
+ /*ATTRS*/
+ );
+ FUNCTOR2D(GridCoGridCoGeom,FrictPhys);
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_GridCoGridCoGeom_FrictPhys_CundallStrack);
+
//!################## Bounds #####################
=== modified file 'pkg/common/InsertionSortCollider.cpp'
--- pkg/common/InsertionSortCollider.cpp 2013-03-06 17:30:45 +0000
+++ pkg/common/InsertionSortCollider.cpp 2013-05-30 20:20:59 +0000
@@ -22,7 +22,6 @@
void InsertionSortCollider::handleBoundInversion(Body::id_t id1, Body::id_t id2, InteractionContainer* interactions, Scene*){
assert(!periodic);
assert(id1!=id2);
- ///fast
if (spatialOverlap(id1,id2) && Collider::mayCollide(Body::byId(id1,scene).get(),Body::byId(id2,scene).get()) && !interactions->found(id1,id2))
interactions->insert(shared_ptr<Interaction>(new Interaction(id1,id2)));
}
=== modified file 'pkg/common/OpenGLRenderer.cpp'
--- pkg/common/OpenGLRenderer.cpp 2013-03-07 18:28:01 +0000
+++ pkg/common/OpenGLRenderer.cpp 2013-05-30 20:13:27 +0000
@@ -113,7 +113,7 @@
// draw periodic cell, if active
void OpenGLRenderer::drawPeriodicCell(){
if(!scene->isPeriodic) return;
- glColor3v(Vector3r(1,1,0));
+ glColor3v(cellColor);
glPushMatrix();
// Vector3r size=scene->cell->getSize();
const Matrix3r& hSize=scene->cell->hSize;
=== modified file 'pkg/common/OpenGLRenderer.hpp'
--- pkg/common/OpenGLRenderer.hpp 2012-04-24 19:42:08 +0000
+++ pkg/common/OpenGLRenderer.hpp 2013-05-30 20:13:27 +0000
@@ -93,6 +93,7 @@
((Vector3r,light2Pos,Vector3r(-130,75,30),,"Position of secondary OpenGL light source in the scene."))
((Vector3r,lightColor,Vector3r(0.6,0.6,0.6),,"Per-color intensity of primary light (RGB)."))
((Vector3r,light2Color,Vector3r(0.5,0.5,0.1),,"Per-color intensity of secondary light (RGB)."))
+ ((Vector3r,cellColor,Vector3r(1,1,0),,"Color of the periodic cell (RGB)."))
((Vector3r,bgColor,Vector3r(.2,.2,.2),,"Color of the background canvas (RGB)"))
((bool,wire,false,,"Render all bodies with wire only (faster)"))
((bool,light1,true,,"Turn light 1 on."))
=== added file 'pkg/dem/BubbleMat.cpp'
--- pkg/dem/BubbleMat.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/BubbleMat.cpp 2013-06-15 19:17:48 +0000
@@ -0,0 +1,78 @@
+#include "BubbleMat.hpp"
+
+
+
+YADE_PLUGIN((BubbleMat)(Ip2_BubbleMat_BubbleMat_BubblePhys)(BubblePhys)(Law2_ScGeom_BubblePhys_Bubble));
+
+
+/********************** Ip2_BubbleMat_BubbleMat_BubblePhys ****************************/
+CREATE_LOGGER(Ip2_BubbleMat_BubbleMat_BubblePhys);
+void Ip2_BubbleMat_BubbleMat_BubblePhys::go(const shared_ptr<Material>& m1, const shared_ptr<Material>& m2, const shared_ptr<Interaction>& interaction){
+ // phys already exists
+ if (interaction->phys) return;
+
+ shared_ptr<BubblePhys> phys(new BubblePhys());
+ interaction->phys = phys;
+ BubbleMat* mat1 = YADE_CAST<BubbleMat*>(m1.get());
+ BubbleMat* mat2 = YADE_CAST<BubbleMat*>(m2.get());
+
+ // averaging over both materials
+ phys->surfaceTension = .5*(mat1->surfaceTension + mat2->surfaceTension);
+}
+
+
+/********************** BubblePhys ****************************/
+CREATE_LOGGER(BubblePhys);
+Real BubblePhys::computeForce(Real penetrationDepth, Real surfaceTension, Real rAvg, int newtonIter, Real newtonTol) {
+ if (penetrationDepth <= 0.0) { return 0.0; }
+
+ Real f,df,ll,retOld,residual;
+ Real c1 = 1./(4*Mathr::PI*surfaceTension);
+ Real c2 = 1./(8*Mathr::PI*surfaceTension*rAvg);
+ Real ret=1./c2;
+ int i = 0;
+ do {
+ retOld = ret;
+ ll = log( ret*c2 );
+ f = penetrationDepth - ret*c1*ll;
+ df = -c1*(ll + 1);
+ ret -= f/df;
+ residual = fabs(ret - retOld)/retOld;
+ if (i++ > newtonIter) {
+ throw runtime_error("BubblePhys::computeForce: no convergence\n");
+ }
+ } while (residual > newtonTol);
+ return ret;
+}
+
+
+
+
+/********************** Law2_ScGeom_BubblePhys_Bubble ****************************/
+CREATE_LOGGER(Law2_ScGeom_BubblePhys_Bubble);
+
+void Law2_ScGeom_BubblePhys_Bubble::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I){
+ ScGeom* geom=static_cast<ScGeom*>(_geom.get());
+ BubblePhys* phys=static_cast<BubblePhys*>(_phys.get());
+
+ if (I->isFresh(scene)) {
+ phys->rAvg = .5*(geom->refR1 + geom->refR2);
+ }
+
+ Real &fN = phys->fN;
+ fN = phys->computeForce(geom->penetrationDepth, phys->surfaceTension, phys->rAvg, phys->newtonIter, phys->newtonTol);
+ Vector3r &normalForce = phys->normalForce;
+ normalForce = fN*geom->normal;
+
+ int id1 = I->getId1(), id2 = I->getId2();
+ if (!scene->isPeriodic) {
+ State* b1 = Body::byId(id1,scene)->state.get();
+ State* b2 = Body::byId(id2,scene)->state.get();
+ applyForceAtContactPoint(-normalForce, geom->contactPoint , id1, b1->se3.position, id2, b2->se3.position);
+ } else {
+ scene->forces.addForce(id1,-normalForce);
+ scene->forces.addForce(id2,normalForce);
+ scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(normalForce));
+ scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(normalForce));
+ }
+}
=== added file 'pkg/dem/BubbleMat.hpp'
--- pkg/dem/BubbleMat.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/BubbleMat.hpp 2013-06-15 19:17:48 +0000
@@ -0,0 +1,75 @@
+#pragma once
+
+#include<yade/core/Material.hpp>
+#include<yade/core/IPhys.hpp>
+#include<yade/pkg/dem/ScGeom.hpp>
+#include<yade/pkg/common/Dispatching.hpp>
+
+namespace py=boost::python;
+
+
+/********************** BubbleMat ****************************/
+class BubbleMat : public Material {
+ public:
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(BubbleMat,Material,"TODO DOC",
+ ((Real,surfaceTension,1/*TODO some realistic value*/,,"TODO DOC"))
+ ,
+ createIndex();
+ density=1; // TODO density default value
+ );
+ REGISTER_CLASS_INDEX(BubbleMat,Material);
+};
+REGISTER_SERIALIZABLE(BubbleMat);
+
+
+/********************** BubblePhys ****************************/
+class BubblePhys : public IPhys {
+ public:
+
+ static Real computeForce(Real penetrationDepth, Real surfaceTension, Real rAvg, int newtonIter, Real newtonTol);
+
+ virtual ~BubblePhys(){};
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(BubblePhys,IPhys,"TODO DOC",
+ ((Vector3r,normalForce,Vector3r::Zero(),,"TODO DOC"))
+ ((Real,surfaceTension,NaN,,"TODO DOC"))
+ ((Real,fN,NaN,,"TODO DOC"))
+ ((Real,rAvg,NaN,,"TODO DOC"))
+ ((Real,newtonIter,50,,"TODO DOC"))
+ ((Real,newtonTol,1e-6,,"TODO DOC"))
+ ,
+ createIndex();
+ ,
+ .def("computeForce",&BubblePhys::computeForce,"TODO DOC")
+ .staticmethod("computeForce")
+ );
+ DECLARE_LOGGER;
+ REGISTER_CLASS_INDEX(BubblePhys,IPhys);
+};
+REGISTER_SERIALIZABLE(BubblePhys);
+
+
+
+/********************** Ip2_BubbleMat_BubbleMat_BubblePhys ****************************/
+class Ip2_BubbleMat_BubbleMat_BubblePhys : public IPhysFunctor{
+ public:
+ virtual void go(const shared_ptr<Material>& m1, const shared_ptr<Material>& m2, const shared_ptr<Interaction>& interaction);
+ FUNCTOR2D(BubbleMat,BubbleMat);
+ DECLARE_LOGGER;
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_BubbleMat_BubbleMat_BubblePhys,IPhysFunctor,"TODO DOC",
+ );
+};
+REGISTER_SERIALIZABLE(Ip2_BubbleMat_BubbleMat_BubblePhys);
+
+
+/********************** Law2_ScGeom_BubblePhys_Bubble ****************************/
+class Law2_ScGeom_BubblePhys_Bubble : public LawFunctor{
+ public:
+ void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* interaction);
+ FUNCTOR2D(GenericSpheresContact,BubblePhys);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_BubblePhys_Bubble,LawFunctor,"TODO DOC",
+ ,
+ /*ctor*/,
+ );
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom_BubblePhys_Bubble);
=== modified file 'pkg/dem/CapillaryPhys.hpp'
--- pkg/dem/CapillaryPhys.hpp 2012-03-27 22:07:35 +0000
+++ pkg/dem/CapillaryPhys.hpp 2013-05-30 20:17:45 +0000
@@ -15,16 +15,21 @@
virtual ~CapillaryPhys();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(CapillaryPhys,FrictPhys,"Physics (of interaction) for Law2_ScGeom_CapillaryPhys_Capillarity.",
+ YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(CapillaryPhys,FrictPhys,"Physics (of interaction) for Law2_ScGeom_CapillaryPhys_Capillarity.",
((bool,meniscus,false,,"Presence of a meniscus if true"))
((bool,isBroken,false,,"If true, capillary force is zero and liquid bridge is inactive."))
- ((Real,CapillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
- ((Real,Vmeniscus,0.,,"Volume of the menicus"))
+ ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
+ ((Real,vMeniscus,0.,,"Volume of the menicus"))
((Real,Delta1,0.,,"Defines the surface area wetted by the meniscus on the smallest grains of radius R1 (R1<R2)"))
((Real,Delta2,0.,,"Defines the surface area wetted by the meniscus on the biggest grains of radius R2 (R1<R2)"))
- ((Vector3r,Fcap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
+ ((Vector3r,fCap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
((short int,fusionNumber,0.,,"Indicates the number of meniscii that overlap with this one"))
- ,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
+ ,/*deprec*/
+ ((Fcap,fCap,"naming convention"))
+ ((CapillaryPressure,capillaryPressure,"naming convention"))
+ ,,
+ createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
+ ,
);
REGISTER_CLASS_INDEX(CapillaryPhys,FrictPhys);
};
=== modified file 'pkg/dem/CapillaryStressRecorder.cpp'
--- pkg/dem/CapillaryStressRecorder.cpp 2011-12-01 16:51:41 +0000
+++ pkg/dem/CapillaryStressRecorder.cpp 2013-05-30 20:17:45 +0000
@@ -46,7 +46,7 @@
Real f1_cap_x=0, f1_cap_y=0, f1_cap_z=0, x1=0, y1=0, z1=0, x2=0, y2=0, z2=0;
Real sig11_cap=0, sig22_cap=0, sig33_cap=0, sig12_cap=0, sig13_cap=0,
- sig23_cap=0, Vwater = 0, CapillaryPressure = 0;
+ sig23_cap=0, Vwater = 0, capillaryPressure = 0;
int j = 0;
InteractionContainer::iterator ii = scene->interactions->begin();
@@ -67,7 +67,7 @@
unsigned int id1 = interaction -> getId1();
unsigned int id2 = interaction -> getId2();
- Vector3r fcap = meniscusParameters->Fcap;
+ Vector3r fcap = meniscusParameters->fCap;
f1_cap_x=fcap[0];
f1_cap_y=fcap[1];
@@ -94,8 +94,8 @@
/// Liquid volume
- Vwater += meniscusParameters->Vmeniscus;
- CapillaryPressure=meniscusParameters->CapillaryPressure;
+ Vwater += meniscusParameters->vMeniscus;
+ capillaryPressure=meniscusParameters->capillaryPressure;
}
@@ -148,7 +148,7 @@
<< lexical_cast<string>(SIG_12_cap) << " "
<< lexical_cast<string>(SIG_13_cap)<< " "
<< lexical_cast<string>(SIG_23_cap)<< " "
- << lexical_cast<string>(CapillaryPressure) << " "
+ << lexical_cast<string>(capillaryPressure) << " "
<< lexical_cast<string>(Sr)<< " "
<< lexical_cast<string>(w)<< " "
<< endl;
=== modified file 'pkg/dem/CapillaryTriaxialTest.cpp'
--- pkg/dem/CapillaryTriaxialTest.cpp 2013-01-07 15:19:31 +0000
+++ pkg/dem/CapillaryTriaxialTest.cpp 2013-05-30 20:17:45 +0000
@@ -338,7 +338,7 @@
// capillary
shared_ptr<Law2_ScGeom_CapillaryPhys_Capillarity> capillaryCohesiveLaw(new Law2_ScGeom_CapillaryPhys_Capillarity);
- capillaryCohesiveLaw->CapillaryPressure = CapillaryPressure;
+ capillaryCohesiveLaw->capillaryPressure = capillaryPressure;
// capillaryCohesiveLaw->fusionDetection = fusionDetection;
// capillaryCohesiveLaw->binaryFusion = binaryFusion;
=== modified file 'pkg/dem/CapillaryTriaxialTest.hpp'
--- pkg/dem/CapillaryTriaxialTest.hpp 2011-01-21 08:14:28 +0000
+++ pkg/dem/CapillaryTriaxialTest.hpp 2013-05-30 20:17:45 +0000
@@ -71,7 +71,7 @@
((string,importFilename,"",,"File with positions and sizes of spheres."))
((string,Key,"",,"A code that is added to output filenames."))
((string,fixedBoxDims,"",,"string that contains some subset (max. 2) of {'x','y','z'} ; contains axes will have box dimension hardcoded, even if box is scaled as mean_radius is prescribed: scaling will be applied on the rest."))
- ((Real,CapillaryPressure,0,,"Define succion in the packing [Pa]. This is the value used in the capillary model."))
+ ((Real,capillaryPressure,0,,"Define succion in the packing [Pa]. This is the value used in the capillary model."))
((bool,water,true,,"activate capillary model"))
((bool,fusionDetection,false,,"test overlaps between liquid bridges on modify forces if overlaps exist"))
((bool,binaryFusion,true,,"Defines how overlapping bridges affect the capillary forces (see :yref:`CapillaryTriaxialTest::fusionDetection`). If binary=true, the force is null as soon as there is an overlap detected, if not, the force is divided by the number of overlaps."))
=== modified file 'pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp 2012-03-27 22:07:35 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp 2013-05-30 20:17:45 +0000
@@ -19,16 +19,20 @@
virtual ~MindlinCapillaryPhys();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(MindlinCapillaryPhys,MindlinPhys,"Adds capillary physics to Mindlin's interaction physics.",
+ YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(MindlinCapillaryPhys,MindlinPhys,"Adds capillary physics to Mindlin's interaction physics.",
((bool,meniscus,false,,"Presence of a meniscus if true"))
((bool,isBroken,false,,"If true, capillary force is zero and liquid bridge is inactive."))
- ((Real,CapillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
- ((Real,Vmeniscus,0.,,"Volume of the menicus"))
+ ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
+ ((Real,vMeniscus,0.,,"Volume of the menicus"))
((Real,Delta1,0.,,"Defines the surface area wetted by the meniscus on the smallest grains of radius R1 (R1<R2)"))
((Real,Delta2,0.,,"Defines the surface area wetted by the meniscus on the biggest grains of radius R2 (R1<R2)"))
- ((Vector3r,Fcap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
+ ((Vector3r,fCap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
((short int,fusionNumber,0.,,"Indicates the number of meniscii that overlap with this one"))
- ,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
+ ,/*deprec*/
+ ((Fcap,fCap,"naming convention"))
+ ((CapillaryPressure,capillaryPressure,"naming convention"))
+ ,,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
+ ,
);
REGISTER_CLASS_INDEX(MindlinCapillaryPhys,MindlinPhys);
};
=== modified file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp 2013-01-29 16:16:15 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp 2013-05-30 20:17:45 +0000
@@ -67,6 +67,7 @@
void Law2_ScGeom_CapillaryPhys_Capillarity::action()
{
if (!scene) cerr << "scene not defined!";
+ if (!capillary) postLoad(*this);//when the script does not define arguments, postLoad is never called
shared_ptr<BodyContainer>& bodies = scene->bodies;
if (fusionDetection && !bodiesMenisciiList.initialized) bodiesMenisciiList.prepare(scene);
@@ -119,7 +120,7 @@
Real D = alpha*((b2->state->pos-b1->state->pos).norm()-(currentContactGeometry->radius1+ currentContactGeometry->radius2)); // scGeom->penetrationDepth could probably be used here?
if ((currentContactGeometry->penetrationDepth>=0)|| D<=0 || createDistantMeniscii) { //||(scene->iter < 1) ) // a simplified way to define meniscii everywhere
- D=0; // defines Fcap when spheres interpenetrate. D<0 leads to wrong interpolation has D<0 has no solution in the interpolation : this is not physically interpretable!! even if, interpenetration << grain radius.
+ D=0; // defines fCap when spheres interpenetrate. D<0 leads to wrong interpolation has D<0 has no solution in the interpolation : this is not physically interpretable!! even if, interpenetration << grain radius.
if (!hertzOn) {
if (fusionDetection && !cundallContactPhysics->meniscus) bodiesMenisciiList.insert((*ii));
cundallContactPhysics->meniscus=true;
@@ -132,10 +133,10 @@
/// Suction (Capillary pressure):
Real Pinterpol = 0;
- if (!hertzOn) Pinterpol = cundallContactPhysics->isBroken ? 0 : CapillaryPressure*(R2/liquidTension);
- else Pinterpol = mindlinContactPhysics->isBroken ? 0 : CapillaryPressure*(R2/liquidTension);
- if (!hertzOn) cundallContactPhysics->CapillaryPressure = CapillaryPressure;
- else mindlinContactPhysics->CapillaryPressure = CapillaryPressure;
+ if (!hertzOn) Pinterpol = cundallContactPhysics->isBroken ? 0 : capillaryPressure*(R2/liquidTension);
+ else Pinterpol = mindlinContactPhysics->isBroken ? 0 : capillaryPressure*(R2/liquidTension);
+ if (!hertzOn) cundallContactPhysics->capillaryPressure = capillaryPressure;
+ else mindlinContactPhysics->capillaryPressure = capillaryPressure;
/// Capillary solution finder:
if ((Pinterpol>=0) && (hertzOn? mindlinContactPhysics->meniscus : cundallContactPhysics->meniscus)) {
@@ -145,17 +146,17 @@
solution(Pinterpol? capillary->Interpolate(R1,R2,Dinterpol, Pinterpol, currentIndexes) : MeniscusParameters());
/// capillary adhesion force
Real Finterpol = solution.F;
- Vector3r Fcap = Finterpol*(2*Mathr::PI*(R2/alpha)*liquidTension)*currentContactGeometry->normal;
- if (!hertzOn) cundallContactPhysics->Fcap = Fcap;
- else mindlinContactPhysics->Fcap = Fcap;
+ Vector3r fCap = Finterpol*(2*Mathr::PI*(R2/alpha)*liquidTension)*currentContactGeometry->normal;
+ if (!hertzOn) cundallContactPhysics->fCap = fCap;
+ else mindlinContactPhysics->fCap = fCap;
/// meniscus volume
Real Vinterpol = solution.V * pow(R2/alpha,3);
if (!hertzOn) {
- cundallContactPhysics->Vmeniscus = Vinterpol;
+ cundallContactPhysics->vMeniscus = Vinterpol;
if (Vinterpol != 0) cundallContactPhysics->meniscus = true;
else cundallContactPhysics->meniscus = false;
} else {
- mindlinContactPhysics->Vmeniscus = Vinterpol;
+ mindlinContactPhysics->vMeniscus = Vinterpol;
if (Vinterpol != 0) mindlinContactPhysics->meniscus = true;
else mindlinContactPhysics->meniscus = false;
}
@@ -190,15 +191,15 @@
short int& fusionNumber = hertzOn?mindlinContactPhysics->fusionNumber:cundallContactPhysics->fusionNumber;
if (binaryFusion) {
if (fusionNumber!=0) { //cerr << "fusion" << endl;
- hertzOn?mindlinContactPhysics->Fcap:cundallContactPhysics->Fcap = Vector3r::Zero();
+ hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap = Vector3r::Zero();
continue;
}
}
//LINEAR VERSION : capillary force is divided by (fusionNumber + 1) - NOTE : any decreasing function of fusionNumber can be considered in fact
- else if (fusionNumber !=0) hertzOn?mindlinContactPhysics->Fcap:cundallContactPhysics->Fcap /= (fusionNumber+1.);
+ else if (fusionNumber !=0) hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap /= (fusionNumber+1.);
}
- scene->forces.addForce((*ii)->getId1(), hertzOn?mindlinContactPhysics->Fcap:cundallContactPhysics->Fcap);
- scene->forces.addForce((*ii)->getId2(),-(hertzOn?mindlinContactPhysics->Fcap:cundallContactPhysics->Fcap));
+ scene->forces.addForce((*ii)->getId1(), hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap);
+ scene->forces.addForce((*ii)->getId2(),-(hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap));
}
}
}
@@ -358,7 +359,7 @@
for (int i=0; i<n_D; i++)
full_data.push_back(TableauD(file));
file.close();
- //cerr << *this; // exemple d'utilisation de la fonction d'ecriture (this est le pointeur vers l'objet courant)
+
}
Tableau::~Tableau()
=== modified file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp 2013-01-29 16:16:15 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp 2013-05-30 20:17:45 +0000
@@ -13,7 +13,6 @@
#include <yade/core/GlobalEngine.hpp>
#include <set>
#include <boost/tuple/tuple.hpp>
-
#include <vector>
#include <list>
#include <utility>
@@ -94,12 +93,15 @@
void action();
void postLoad(Law2_ScGeom_CapillaryPhys_Capillarity&);
- YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CapillaryPhys_Capillarity,GlobalEngine,"This law allows one to take into account capillary forces/effects between spheres coming from the presence of interparticular liquid bridges (menisci).\n\nThe control parameter is the capillary pressure (or suction) Uc = ugas - Uliquid. Liquid bridges properties (volume V, extent over interacting grains delta1 and delta2) are computed as a result of the defined capillary pressure and of the interacting geometry (spheres radii and interparticular distance).\n\nReferences: in english [Scholtes2009b]_; more detailed, but in french [Scholtes2009d]_.\n\nThe law needs ascii files M(r=i) with i=R1/R2 to work (see https://yade-dem.org/index.php/CapillaryTriaxialTest). These ASCII files contain a set of results from the resolution of the Laplace-Young equation for different configurations of the interacting geometry.\n\nIn order to allow capillary forces between distant spheres, it is necessary to enlarge the bounding boxes using :yref:`Bo1_Sphere_Aabb::aabbEnlargeFactor` and make the Ig2 define define distant interactions via :yref:`interactionDetectionFactor<Ig2_Sphere_Sphere_ScGeom::interactionDetectionFactor>`. It is also necessary to disable interactions removal by the constitutive law (:yref:`Law2<Law2_ScGeom_FrictPhys_CundallStrack::neverErase>=True`). The only combinations of laws supported are currently capillary law + :yref:`Law2_ScGeom_FrictPhys_CundallStrack` and capillary law + :yref:`Law2_ScGeom_MindlinPhys_Mindlin` (and the other variants of Hertz-Mindlin).\n\nSee CapillaryPhys-example.py for an example script.",
- ((Real,CapillaryPressure,0.,,"Value of the capillary pressure Uc defines as Uc=Ugas-Uliquid"))
+ YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(Law2_ScGeom_CapillaryPhys_Capillarity,GlobalEngine,"This law allows one to take into account capillary forces/effects between spheres coming from the presence of interparticular liquid bridges (menisci).\n\nThe control parameter is the capillary pressure (or suction) Uc = ugas - Uliquid. Liquid bridges properties (volume V, extent over interacting grains delta1 and delta2) are computed as a result of the defined capillary pressure and of the interacting geometry (spheres radii and interparticular distance).\n\nReferences: in english [Scholtes2009b]_; more detailed, but in french [Scholtes2009d]_.\n\nThe law needs ascii files M(r=i) with i=R1/R2 to work (see https://yade-dem.org/index.php/CapillaryTriaxialTest). These ASCII files contain a set of results from the resolution of the Laplace-Young equation for different configurations of the interacting geometry.\n\nIn order to allow capillary forces between distant spheres, it is necessary to enlarge the bounding boxes using :yref:`Bo1_Sphere_Aabb::aabbEnlargeFactor` and make the Ig2 define define distant interactions via :yref:`interactionDetectionFactor<Ig2_Sphere_Sphere_ScGeom::interactionDetectionFactor>`. It is also necessary to disable interactions removal by the constitutive law (:yref:`Law2<Law2_ScGeom_FrictPhys_CundallStrack::neverErase>=True`). The only combinations of laws supported are currently capillary law + :yref:`Law2_ScGeom_FrictPhys_CundallStrack` and capillary law + :yref:`Law2_ScGeom_MindlinPhys_Mindlin` (and the other variants of Hertz-Mindlin).\n\nSee CapillaryPhys-example.py for an example script.",
+ ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Uc=Ugas-Uliquid"))
((bool,fusionDetection,false,,"If true potential menisci overlaps are checked"))
((bool,binaryFusion,true,,"If true, capillary forces are set to zero as soon as, at least, 1 overlap (menisci fusion) is detected"))
((bool,hertzOn,false,,"|yupdate| true if hertz model is used"))
- ((bool,createDistantMeniscii,false,,"Generate meniscii between distant spheres? Else only maintain the existing one. For modeling a wetting path this flag should always be false. For a drying path it should be true for one step (initialization) then false, as in the logic of [Scholtes2009c]_"))
+ ((bool,createDistantMeniscii,false,,"Generate meniscii between distant spheres? Else only maintain the existing ones. For modeling a wetting path this flag should always be false. For a drying path it should be true for one step (initialization) then false, as in the logic of [Scholtes2009c]_"))
+ ,/*deprec*/
+ ((CapillaryPressure,capillaryPressure,"naming convention"))
+ ,,,
);
};
=== modified file 'pkg/dem/SampleCapillaryPressureEngine.cpp'
--- pkg/dem/SampleCapillaryPressureEngine.cpp 2013-02-07 18:12:42 +0000
+++ pkg/dem/SampleCapillaryPressureEngine.cpp 2013-05-30 20:17:45 +0000
@@ -60,12 +60,12 @@
if (scene->iter % 100 == 0) cerr << "pressure variation!!" << endl;
if ((Pressure>=0) && (Pressure<=1000000000)) Pressure += PressureVariation;
- capillaryCohesiveLaw->CapillaryPressure = Pressure;
+ capillaryCohesiveLaw->capillaryPressure = Pressure;
capillaryCohesiveLaw->fusionDetection = fusionDetection;
capillaryCohesiveLaw->binaryFusion = binaryFusion;
}
- else { capillaryCohesiveLaw->CapillaryPressure = Pressure;
+ else { capillaryCohesiveLaw->capillaryPressure = Pressure;
capillaryCohesiveLaw->fusionDetection = fusionDetection;
capillaryCohesiveLaw->binaryFusion = binaryFusion;}
if (scene->iter % 100 == 0) cerr << "capillary pressure = " << Pressure << endl;
=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp 2013-03-19 09:58:00 +0000
+++ pkg/dem/Shop.cpp 2013-05-30 20:19:43 +0000
@@ -19,6 +19,7 @@
#include<yade/pkg/common/Sphere.hpp>
#include<yade/pkg/common/ElastMat.hpp>
#include<yade/pkg/dem/ViscoelasticPM.hpp>
+#include<yade/pkg/dem/CapillaryPhys.hpp>
#include<yade/pkg/common/Bo1_Sphere_Aabb.hpp>
#include<yade/pkg/common/Bo1_Box_Aabb.hpp>
@@ -778,6 +779,23 @@
return stressTensor/volume;
}
+Matrix3r Shop::getCapillaryStress(Real volume){
+ Scene* scene=Omega::instance().getScene().get();
+ if (volume==0) volume = scene->isPeriodic?scene->cell->hSize.determinant():1;
+ Matrix3r stressTensor = Matrix3r::Zero();
+ const bool isPeriodic = scene->isPeriodic;
+ FOREACH(const shared_ptr<Interaction>&I, *scene->interactions){
+ if (!I->isReal()) continue;
+ shared_ptr<Body> b1 = Body::byId(I->getId1(),scene);
+ shared_ptr<Body> b2 = Body::byId(I->getId2(),scene);
+ CapillaryPhys* nsi=YADE_CAST<CapillaryPhys*> ( I->phys.get() );
+ Vector3r branch=b1->state->pos -b2->state->pos;
+ if (isPeriodic) branch-= scene->cell->hSize*I->cellDist.cast<Real>();
+ stressTensor += (nsi->fCap)*branch.transpose();
+ }
+ return stressTensor/volume;
+}
+
/*
Matrix3r Shop::stressTensorOfPeriodicCell(bool smallStrains){
Scene* scene=Omega::instance().getScene().get();
=== modified file 'pkg/dem/Shop.hpp'
--- pkg/dem/Shop.hpp 2013-01-07 15:19:31 +0000
+++ pkg/dem/Shop.hpp 2013-05-30 20:19:43 +0000
@@ -121,6 +121,7 @@
//! Function to compute overall ("macroscopic") stress.
static Matrix3r getStress(Real volume=0);
+ static Matrix3r getCapillaryStress(Real volume=0);
static Matrix3r stressTensorOfPeriodicCell() { LOG_WARN("Shop::stressTensorOfPeriodicCelli is DEPRECATED: use getStress instead"); return Shop::getStress(); }
//! This version is restricted to periodic BCs and Dem3Dof
static Matrix3r stressTensorOfPeriodicCell(bool smallStrains=true);
=== added file 'pkg/dem/ViscoelasticCapillarPM.cpp'
--- pkg/dem/ViscoelasticCapillarPM.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/ViscoelasticCapillarPM.cpp 2013-06-12 14:35:50 +0000
@@ -0,0 +1,138 @@
+#include"ViscoelasticPM.hpp"
+#include<yade/core/State.hpp>
+#include<yade/pkg/dem/ScGeom.hpp>
+#include<yade/core/Omega.hpp>
+#include<yade/core/Scene.hpp>
+#include<yade/pkg/common/Sphere.hpp>
+
+Real Law2_ScGeom_ViscElPhys_Basic::calculateCapillarForce(const ScGeom& geom, ViscElPhys& phys) {
+ Real fC = 0.0;
+
+ /* Capillar
+ * Some equations have constants, which can be calculated only once per contact.
+ * No need to recalculate them at each step.
+ * It needs to be fixed.
+ *
+ */
+
+ if (phys.CapillarType == "Weigert") {
+ /* Capillar model from [Weigert1999]
+ */
+ Real R = phys.R;
+ Real a = -geom.penetrationDepth;
+ Real Ca = (1.0 + 6.0*a/(R*2.0)); // [Weigert1999], equation (16)
+ Real Ct = (1.0 + 1.1*sin(phys.theta)); // [Weigert1999], equation (17)
+
+ /*
+ Real Eps = 0.36; // Porosity
+ Real fi = phys.Vb/(2.0*M_PI/6.0*pow(R*2.0,3.)); // [Weigert1999], equation (13)
+ Real S = M_PI*(1-Eps)/(pow(Eps, 2.0))*fi; // [Weigert1999], equation (14)
+ Real beta = asin(pow(((S/0.36)*(pow(Eps, 2.0)/(1-Eps))*(1.0/Ca)*(1.0/Ct)), 1.0/4.0)); // [Weigert1999], equation (19)
+ */
+
+ Real beta = asin(pow(phys.Vb/(0.12*Ca*Ct*pow(2.0*R, 3.0)), 1.0/4.0)); // [Weigert1999], equation (15), against Vb
+
+ Real r1 = (2.0*R*(1-cos(beta)) + a)/(2.0*cos(beta+phys.theta)); // [Weigert1999], equation (5)
+ Real r2 = R*sin(beta) + r1*(sin(beta+phys.theta)-1); // [Weigert1999], equation (6)
+ Real Pk = phys.gamma*(1/r1 - 1/r2); /* [Weigert1999], equation (22),
+ * see also a sentence over the equation
+ * "R1 was taken as positive and R2 was taken as negative"
+ */
+
+ //fC = M_PI*2.0*R*phys.gamma/(1+tan(0.5*beta)); // [Weigert1999], equation (23), [Fisher]
+
+ fC = M_PI/4.0*pow((2.0*R),2.0)*pow(sin(beta),2.0)*Pk +
+ phys.gamma*M_PI*2.0*R*sin(beta)*sin(beta+phys.theta); // [Weigert1999], equation (21)
+
+ } else if (phys.CapillarType == "Willett_numeric") {
+
+ /* Capillar model from [Willett2000]
+ */
+
+ Real R = phys.R;
+ Real s = -geom.penetrationDepth/2.0;
+ Real Vb = phys.Vb;
+
+ Real VbS = Vb/(R*R*R);
+ Real Th1 = phys.theta;
+ Real Th2 = phys.theta*phys.theta;
+ Real Gamma = phys.gamma;
+
+ /*
+ * [Willett2000], equations in Anhang
+ */
+ Real f1 = (-0.44507 + 0.050832*Th1 - 1.1466*Th2) +
+ (-0.1119 - 0.000411*Th1 - 0.1490*Th2) * log(VbS) +
+ (-0.012101 - 0.0036456*Th1 - 0.01255*Th2) *log(VbS)*log(VbS) +
+ (-0.0005 - 0.0003505*Th1 - 0.00029076*Th2) *log(VbS)*log(VbS)*log(VbS);
+
+ Real f2 = (1.9222 - 0.57473*Th1 - 1.2918*Th2) +
+ (-0.0668 - 0.1201*Th1 - 0.22574*Th2) * log(VbS) +
+ (-0.0013375 - 0.0068988*Th1 - 0.01137*Th2) *log(VbS)*log(VbS);
+
+ Real f3 = (1.268 - 0.01396*Th1 - 0.23566*Th2) +
+ (0.198 + 0.092*Th1 - 0.06418*Th2) * log(VbS) +
+ (0.02232 + 0.02238*Th1 - 0.009853*Th2) *log(VbS)*log(VbS) +
+ (0.0008585 + 0.001318*Th1 - 0.00053*Th2) *log(VbS)*log(VbS)*log(VbS);
+
+ Real f4 = (-0.010703 + 0.073776*Th1 - 0.34742*Th2) +
+ (0.03345 + 0.04543*Th1 - 0.09056*Th2) * log(VbS) +
+ (0.0018574 + 0.004456*Th1 - 0.006257*Th2) *log(VbS)*log(VbS);
+
+ Real sPl = s/sqrt(Vb/R);
+
+ Real lnFS = f1 - f2*exp(f3*log(sPl) + f4*log(sPl)*log(sPl));
+ Real FS = exp(lnFS);
+
+ fC = FS * 2.0 * M_PI* R * Gamma;
+ } else if (phys.CapillarType == "Willett_analytic") {
+ /* Capillar model from Willet [Willett2000] (analytical solution), but
+ * used also in the work of Herminghaus [Herminghaus2005]
+ */
+
+ Real R = phys.R;
+ Real Gamma = phys.gamma;
+ Real s = -geom.penetrationDepth;
+ Real Vb = phys.Vb;
+
+ /*
+
+ Real sPl = s/sqrt(Vb/R); // [Herminghaus2005], equation (sentence between (7) and (8))
+ fC = 2.0 * M_PI* R * Gamma * cos(phys.theta)/(1 + 1.05*sPl + 2.5 *sPl * sPl); // [Herminghaus2005], equation (7)
+
+ */
+
+ Real sPl = (s/2.0)/sqrt(Vb/R); // [Willett2000], equation (sentence after (11)), s - half-separation, so s*2.0
+ Real f_star = cos(phys.theta)/(1 + 2.1*sPl + 10.0 * pow(sPl, 2.0)); // [Willett2000], equation (12)
+ fC = f_star * (2*M_PI*R*Gamma); // [Willett2000], equation (13), against F
+
+ } else if (phys.CapillarType == "Rabinovich") {
+ /* Capillar model from Rabinovich [Rabinov2005]
+ */
+
+ Real R = phys.R;
+ Real Gamma = phys.gamma;
+ Real H = -geom.penetrationDepth;
+ Real V = phys.Vb;
+
+ Real alpha = 0.0;
+ Real dsp = 0.0;
+ if (H!=0.0) {
+ alpha = sqrt(H/R*(-1+ sqrt(1 + 2.0*V/(M_PI*R*H*H)))); // [Rabinov2005], equation (A3)
+
+ dsp = H/2.0*(-1.0 + sqrt(1.0 + 2.0*V/(M_PI*R*H*H))); // [Rabinov2005], equation (20)
+
+ fC = -(2*M_PI*R*Gamma*cos(phys.theta))/(1+(H/(2*dsp))) -
+ 2*M_PI*R*Gamma*sin(alpha)*sin(phys.theta + alpha); // [Rabinov2005], equation (19)
+ } else {
+ fC = -(2*M_PI*R*Gamma*cos(phys.theta)) -
+ 2*M_PI*R*Gamma*sin(alpha)*sin(phys.theta + alpha); // [Rabinov2005], equation (19)
+ }
+
+ fC *=-1;
+
+ } else {
+ throw runtime_error("CapillarType is unknown, please, use only Willett_numeric, Willett_analytic, Weigert or Rabinovich");
+ }
+ return fC;
+}
=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp 2013-05-28 06:22:52 +0000
+++ pkg/dem/ViscoelasticPM.cpp 2013-06-12 14:35:50 +0000
@@ -20,13 +20,20 @@
if(interaction->phys) return;
ViscElMat* mat1 = static_cast<ViscElMat*>(b1.get());
ViscElMat* mat2 = static_cast<ViscElMat*>(b2.get());
- const Real mass1 = Body::byId(interaction->getId1())->state->mass;
- const Real mass2 = Body::byId(interaction->getId2())->state->mass;
+ Real mass1 = 1.0;
+ Real mass2 = 1.0;
+
+ if (mat1->massMultiply and mat2->massMultiply) {
+ mass1 = Body::byId(interaction->getId1())->state->mass;
+ mass2 = Body::byId(interaction->getId2())->state->mass;
+ }
+
const Real kn1 = mat1->kn*mass1; const Real cn1 = mat1->cn*mass1;
const Real ks1 = mat1->ks*mass1; const Real cs1 = mat1->cs*mass1;
+ const Real mR1 = mat1->mR; const Real mR2 = mat2->mR;
+ const int mRtype1 = mat1->mRtype; const int mRtype2 = mat2->mRtype;
const Real kn2 = mat2->kn*mass2; const Real cn2 = mat2->cn*mass2;
const Real ks2 = mat2->ks*mass2; const Real cs2 = mat2->cs*mass2;
-
ViscElPhys* phys = new ViscElPhys();
@@ -41,12 +48,24 @@
phys->ks = 0;
}
+ if ((mR1>0) or (mR2>0)) {
+ phys->mR = 2.0/( ((mR1>0)?1/mR1:0) + ((mR2>0)?1/mR2:0) );
+ } else {
+ phys->mR = 0;
+ }
+
phys->cn = (cn1?1/cn1:0) + (cn2?1/cn2:0); phys->cn = phys->cn?1/phys->cn:0;
phys->cs = (cs1?1/cs1:0) + (cs2?1/cs2:0); phys->cs = phys->cs?1/phys->cs:0;
-
+
phys->tangensOfFrictionAngle = std::tan(std::min(mat1->frictionAngle, mat2->frictionAngle));
phys->shearForce = Vector3r(0,0,0);
+ if ((mRtype1 != mRtype2) or (mRtype1>2) or (mRtype2>2) or (mRtype1<1) or (mRtype2<1) ) {
+ throw runtime_error("mRtype should be equal for both materials and have the values 1 or 2!");
+ } else {
+ phys->mRtype = mRtype1;
+ }
+
if (mat1->Capillar and mat2->Capillar) {
if (mat1->Vb == mat2->Vb) {
phys->Vb = mat1->Vb;
@@ -141,7 +160,7 @@
const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel;
const Real normalVelocity = geom.normal.dot(relativeVelocity);
const Vector3r shearVelocity = relativeVelocity-normalVelocity*geom.normal;
-
+
// As Chiara Modenese suggest, we store the elastic part
// and then add the viscous part if we pass the Mohr-Coulomb criterion.
// See http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg01391.html
@@ -149,8 +168,20 @@
Vector3r shearForceVisc = Vector3r::Zero(); // the viscous shear damping haven't a history because it is a function of the instant velocity
phys.normalForce = ( phys.kn * geom.penetrationDepth + phys.cn * normalVelocity ) * geom.normal;
- //phys.prevNormal = geom.normal;
+
+ Vector3r momentResistance = Vector3r::Zero();
+ if (phys.mR>0.0) {
+ const Vector3r relAngVel = de1.angVel - de2.angVel;
+ relAngVel.normalized();
+
+ if (phys.mRtype == 1) {
+ momentResistance = -phys.mR*phys.normalForce.norm()*relAngVel; // [Zhou1999536], equation (3)
+ } else if (phys.mRtype == 2) {
+ momentResistance = -phys.mR*(c1x.cross(de1.angVel) - c2x.cross(de2.angVel)).norm()*phys.normalForce.norm()*relAngVel; // [Zhou1999536], equation (4)
+ }
+ }
+
const Real maxFs = phys.normalForce.squaredNorm() * std::pow(phys.tangensOfFrictionAngle,2);
if( shearForce.squaredNorm() > maxFs )
{
@@ -167,118 +198,10 @@
}
if (I->isActive) {
- //std::cerr<<"Contact: "<<phys.normalForce<<std::endl;
const Vector3r f = phys.normalForce + shearForce + shearForceVisc;
addForce (id1,-f,scene);
addForce (id2, f,scene);
- addTorque(id1,-c1x.cross(f),scene);
- addTorque(id2, c2x.cross(f),scene);
+ addTorque(id1,-c1x.cross(f)+momentResistance,scene);
+ addTorque(id2, c2x.cross(f)-momentResistance,scene);
}
}
-
-Real Law2_ScGeom_ViscElPhys_Basic::calculateCapillarForce(const ScGeom& geom, ViscElPhys& phys) {
- Real fC = 0.0;
-
- /* Capillar
- * Some equations have constants, which can be calculated only once per contact.
- * No need to recalculate them at each step.
- * It needs to be fixed.
- *
- */
-
- if (phys.CapillarType == "Weigert") {
- /* Capillar model from [Weigert1999]
- */
- Real R = phys.R;
- Real a = -geom.penetrationDepth;
- Real Ca = (1.0 + 6.0*a/(R*2.0)); // [Weigert1999], equation (16)
- Real Ct = (1.0 + 1.1*sin(phys.theta)); // [Weigert1999], equation (17)
-
- /*
- Real Eps = 0.36; // Porosity
- Real fi = phys.Vb/(2.0*M_PI/6.0*pow(R*2.0,3.)); // [Weigert1999], equation (13)
- Real S = M_PI*(1-Eps)/(pow(Eps, 2.0))*fi; // [Weigert1999], equation (14)
- Real beta = asin(pow(((S/0.36)*(pow(Eps, 2.0)/(1-Eps))*(1.0/Ca)*(1.0/Ct)), 1.0/4.0)); // [Weigert1999], equation (19)
- */
-
- Real beta = asin(pow(phys.Vb/(0.12*Ca*Ct*pow(2.0*R, 3.0)), 1.0/4.0)); // [Weigert1999], equation (15), against Vb
-
- Real r1 = (2.0*R*(1-cos(beta)) + a)/(2.0*cos(beta+phys.theta)); // [Weigert1999], equation (5)
- Real r2 = R*sin(beta) + r1*(sin(beta+phys.theta)-1); // [Weigert1999], equation (6)
- Real Pk = phys.gamma*(1/r1 - 1/r2); /* [Weigert1999], equation (22),
- * see also a sentence over the equation
- * "R1 was taken as positive and R2 was taken as negative"
- */
-
- //fC = M_PI*2.0*R*phys.gamma/(1+tan(0.5*beta)); // [Weigert1999], equation (23), [Fisher]
-
- fC = M_PI/4.0*pow((2.0*R),2.0)*pow(sin(beta),2.0)*Pk +
- phys.gamma*M_PI*2.0*R*sin(beta)*sin(beta+phys.theta); // [Weigert1999], equation (21)
-
- } else if (phys.CapillarType == "Willett_numeric") {
-
- /* Capillar model from [Willett2000]
- */
-
- Real R = phys.R;
- Real s = -geom.penetrationDepth/2.0;
- Real Vb = phys.Vb;
-
- Real VbS = Vb/(R*R*R);
- Real Th1 = phys.theta;
- Real Th2 = phys.theta*phys.theta;
- Real Gamma = phys.gamma;
-
- /*
- * [Willett2000], equations in Anhang
- */
- Real f1 = (-0.44507 + 0.050832*Th1 - 1.1466*Th2) +
- (-0.1119 - 0.000411*Th1 - 0.1490*Th2) * log(VbS) +
- (-0.012101 - 0.0036456*Th1 - 0.01255*Th2) *log(VbS)*log(VbS) +
- (-0.0005 - 0.0003505*Th1 - 0.00029076*Th2) *log(VbS)*log(VbS)*log(VbS);
-
- Real f2 = (1.9222 - 0.57473*Th1 - 1.2918*Th2) +
- (-0.0668 - 0.1201*Th1 - 0.22574*Th2) * log(VbS) +
- (-0.0013375 - 0.0068988*Th1 - 0.01137*Th2) *log(VbS)*log(VbS);
-
- Real f3 = (1.268 - 0.01396*Th1 - 0.23566*Th2) +
- (0.198 + 0.092*Th1 - 0.06418*Th2) * log(VbS) +
- (0.02232 + 0.02238*Th1 - 0.009853*Th2) *log(VbS)*log(VbS) +
- (0.0008585 + 0.001318*Th1 - 0.00053*Th2) *log(VbS)*log(VbS)*log(VbS);
-
- Real f4 = (-0.010703 + 0.073776*Th1 - 0.34742*Th2) +
- (0.03345 + 0.04543*Th1 - 0.09056*Th2) * log(VbS) +
- (0.0018574 + 0.004456*Th1 - 0.006257*Th2) *log(VbS)*log(VbS);
-
- Real sPl = s/sqrt(Vb/R);
-
- Real lnFS = f1 - f2*exp(f3*log(sPl) + f4*log(sPl)*log(sPl));
- Real FS = exp(lnFS);
-
- fC = FS * 2.0 * M_PI* R * Gamma;
- } else if (phys.CapillarType == "Willett_analytic") {
- /* Capillar model from Willet [Willett2000] (analytical solution), but
- * used also in the work of Herminghaus [Herminghaus2005]
- */
-
- Real R = phys.R;
- Real Gamma = phys.gamma;
- Real s = -geom.penetrationDepth;
- Real Vb = phys.Vb;
-
- /*
-
- Real sPl = s/sqrt(Vb/R); // [Herminghaus2005], equation (sentence between (7) and (8))
- fC = 2.0 * M_PI* R * Gamma * cos(phys.theta)/(1 + 1.05*sPl + 2.5 *sPl * sPl); // [Herminghaus2005], equation (7)
-
- */
-
- Real sPl = (s/2.0)/sqrt(Vb/R); // [Willett2000], equation (sentence after (11)), s - half-separation, so s*2.0
- Real f_star = cos(phys.theta)/(1 + 2.1*sPl + 10.0 * pow(sPl, 2.0)); // [Willett2000], equation (12)
- fC = f_star * (2*M_PI*R*Gamma); // [Willett2000], equation (13), against F
-
- } else {
- throw runtime_error("CapillarType is unknown, please, use only Willett_numeric, Willett_analytic or Weigert");
- }
- return fC;
-}
=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp 2013-05-28 06:22:52 +0000
+++ pkg/dem/ViscoelasticPM.hpp 2013-06-12 14:35:49 +0000
@@ -22,11 +22,14 @@
((Real,ks,NaN,,"Shear elastic stiffness"))
((Real,cs,NaN,,"Shear viscous constant"))
((Real,frictionAngle,NaN,,"Friction angle [rad]"))
+ ((bool,massMultiply,true,,"Stiffness and viscosity are multiplied by the reduced mass. If massMultiply=false, these parameter are set explicitly without mass multiplication"))
+ ((Real,mR,0.0,,"Rolling resistance, see [Zhou1999536]_."))
+ ((unsigned int,mRtype,1,,"Rolling resistance type, see [Zhou1999536]_. mRtype=1 - equation (3) in [Zhou1999536]_; mRtype=2 - equation (4) in [Zhou1999536]_."))
((bool,Capillar,false,,"True, if capillar forces need to be added."))
((Real,Vb,NaN,,"Liquid bridge volume [m^3]"))
((Real,gamma,NaN,,"Surface tension [N/m]"))
((Real,theta,NaN,,"Contact angle [°]"))
- ((std::string,CapillarType,"",,"Different types of capillar interaction: Willett_numeric, Willett_analytic [Willett2000]_ , Weigert [Weigert1999]_ ")),
+ ((std::string,CapillarType,"",,"Different types of capillar interaction: Willett_numeric, Willett_analytic [Willett2000]_ , Weigert [Weigert1999]_ , Rabinovich [Rabinov2005]_ ")),
createIndex();
);
REGISTER_CLASS_INDEX(ViscElMat,Material);
@@ -41,13 +44,15 @@
YADE_CLASS_BASE_DOC_ATTRS_CTOR(ViscElPhys,FrictPhys,"IPhys created from :yref:`ViscElMat`, for use with :yref:`Law2_ScGeom_ViscElPhys_Basic`.",
((Real,cn,NaN,,"Normal viscous constant"))
((Real,cs,NaN,,"Shear viscous constant"))
+ ((Real,mR,0.0,,"Rolling resistance, see [Zhou1999536]_."))
+ ((unsigned int,mRtype,1,,"Rolling resistance type, see [Zhou1999536]_. mRtype=1 - equation (3) in [Zhou1999536]_; mRtype=2 - equation (4) in [Zhou1999536]_"))
((bool,Capillar,false,,"True, if capillar forces need to be added."))
((bool,liqBridgeCreated,false,,"Whether liquid bridge was created, only after a normal contact of spheres"))
((Real,sCrit,false,,"Critical bridge length [m]"))
((Real,Vb,NaN,,"Liquid bridge volume [m^3]"))
((Real,gamma,NaN,,"Surface tension [N/m]"))
((Real,theta,NaN,,"Contact angle [rad]"))
- ((std::string,CapillarType,"",,"Different types of capillar interaction: Willett, Weigert, Herminghaus")),
+ ((std::string,CapillarType,"",,"Different types of capillar interaction: Willett_numeric, Willett_analytic, Weigert, Rabinovich")),
createIndex();
)
};
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2013-04-02 13:57:36 +0000
+++ py/_utils.cpp 2013-05-31 13:27:13 +0000
@@ -533,6 +533,7 @@
py::def("fabricTensor",Shop__fabricTensor,(py::args("splitTensor")=false,py::args("revertSign")=false,py::args("thresholdForce")=NaN),"Compute the fabric tensor of the periodic cell. The original paper can be found in [Satake1982]_.\n\n:param bool splitTensor: split the fabric tensor into two parts related to the strong and weak contact forces respectively.\n\n:param bool revertSign: it must be set to true if the contact law's convention takes compressive forces as positive.\n\n:param Real thresholdForce: if the fabric tensor is split into two parts, a threshold value can be specified otherwise the mean contact force is considered by default. It is worth to note that this value has a sign and the user needs to set it according to the convention adopted for the contact law. To note that this value could be set to zero if one wanted to make distinction between compressive and tensile forces.");
py::def("bodyStressTensors",Shop__getStressLWForEachBody,"Compute and return a table with per-particle stress tensors. Each tensor represents the average stress in one particle, obtained from the contour integral of applied load as detailed below. This definition is considering each sphere as a continuum. It can be considered exact in the context of spheres at static equilibrium, interacting at contact points with negligible volume changes of the solid phase (this last assumption is not restricting possible deformations and volume changes at the packing scale).\n\nProof:\n\nFirst, we remark the identity: $\\sigma_{ij}=\\delta_{ik}\\sigma_{kj}=x_{i,k}\\sigma_{kj}=(x_{i}\\sigma_{kj})_{,k}-x_{i}\\sigma_{kj,k}$.\n\nAt equilibrium, the divergence of stress is null: $\\sigma_{kj,k}=\\vec{0}$. Consequently, after divergence theorem: $\\frac{1}{V}\\int_V \\sigma_{ij}dV = \\frac{1}{V}\\int_V (x_{i}\\sigma_{kj})_{,k}dV = \\frac{1}{V}\\int_{\\partial V}x_i\\sigma_{kj}n_kdS = \\frac{1}{V}\\sum_bx_i^bf_j^b$.\n\nThe last equality is implicitely based on the representation of external loads as Dirac distributions whose zeros are the so-called *contact points*: 0-sized surfaces on which the *contact forces* are applied, located at $x_i$ in the deformed configuration.\n\nA weighted average of per-body stresses will give the average stress inside the solid phase. There is a simple relation between the stress inside the solid phase and the stress in an equivalent continuum in the absence of fluid pressure. For porosity $n$, the relation reads: $\\sigma_{ij}^{equ.}=(1-n)\\sigma_{ij}^{solid}$.");
py::def("getStress",Shop::getStress,(py::args("volume")=0),"Compute and return Love-Weber stress tensor:\n\n $\\sigma_{ij}=\\frac{1}{V}\\sum_b l_i^b f_j^b$, where the sum is over all interactions, with $l$ the branch vector (joining centers of the bodies) and $f$ is the contact force. $V$ can be passed to the function. If it is not, it will be equal to one in non-periodic cases, or equal to the volume of the cell in periodic cases.");
+ py::def("getCapillaryStress",Shop::getCapillaryStress,(py::args("volume")=0),"Compute and return Love-Weber capillary stress tensor:\n\n $\\sigma^{cap}_{ij}=\\frac{1}{V}\\sum_b l_i^b f^{cap,b}_j$, where the sum is over all interactions, with $l$ the branch vector (joining centers of the bodies) and $f^{cap}$ is the capillary force. $V$ can be passed to the function. If it is not, it will be equal to one in non-periodic cases, or equal to the volume of the cell in periodic cases. Only the CapillaryPhys interaction type is supported presently.");
py::def("getBodyIdsContacts",Shop__getBodyIdsContacts,(py::args("bodyID")=0),"Get a list of body-ids, which contacts the given body.");
py::def("maxOverlapRatio",maxOverlapRatio,"Return maximum overlap ration in interactions (with :yref:`ScGeom`) of two :yref:`spheres<Sphere>`. The ratio is computed as $\\frac{u_N}{2(r_1 r_2)/r_1+r_2}$, where $u_N$ is the current overlap distance and $r_1$, $r_2$ are radii of the two spheres in contact.");
py::def("shiftBodies",shiftBodies,(py::arg("ids"),py::arg("shift")),"Shifts bodies listed in ids without updating their velocities.");
