yade-dev team mailing list archive

[Václav S(milauer <eudoxos@xxxxxxxx>]

>>  I don't use torsion and bending, but for the shear part,
>> I never got the incremental algo to match with what small-strain FEM
>> code computes for strains (production code, therefore assumed to be
>> correct).   
> It is not the first time I read this but I must admit I never
> understodd what it means. What did you compare exactly?
Tension-compression test with 2000 particles of concrete, I gave some
strain to the specimen (about 2e-6) and let it stabilize in DEM. In the
(static) FE code, the same geometrical configuration was used,
interactions are special kind of elements (with the same constitutive
law) and spheres are nodes; I prescribe the same displacement to the
boundary spheres (nodes).

Then I compared strains and stresses on interactions (I picked a few of
them randomly) and they matched (relative difference 1e-5 and such,
perhaps due to slight geometrical non-linearity even at such small
strain). That was for the elastic part. For non-elastic, I compared
strain-stress curve and it nicely matched for DE/FE, except for dynamic
effects that vanished for extremely slow loading in DE.
>> Besides that, implementing geometry algorithms inside constitutive laws
>> is conceptually flawed, for 2 reasons.
>>
>> 1. Constitutive law should operate just on interaction geometry and
>> interaction physics. If it accesses rootBody->bodies, it is a design
>> problem.   
> Agreed : we have a design problem. But I see no obvious solution for
> now (see below).
>>  
> Well... my main point was to save your brain-time providing
> Janek-brained equations.
Eh, I've done it before Christmas. Adding bending/torsion (add
initRelOri12 and methods to compute that) once shear was working was
done in an hour, including testing, I think. I don't need
bending/torsion myself, but I added it to ElasticContactLaw2: it has
normal, shear, bending and torsion parts and has about 20 lines of code
- it is really just demo, but you get the idea how simple it was to
write such law once the geometry things were in place.
> On the c++ part, one solution would be to inherit from
> CohesiveFrictionalContactLaw in some way, or develop the law itself to
> make it do what you need. This law is the most general in Yade
> currently, you can go from the simplest normal-elastic interaction to
> a complex force-moment creep law by playing with the parameters. I'm
> sure it is no big deal to add you own constitutive behaviour in this
> class. Of course, it would be no big deal probably to do the opposite
> : move the CohesiveFrictional stuff to your framework...
You cannot inherit from constitutive laws we have now easily, at least I
don't see how, since it has the loop inside ::action. I _could_ of
course add everything to CohesiveFrictionalContactLaw, but, heck, if I
add viscoplasticity and damage and whatnot, isn't the law getting too
much stuff? To me it seems like writing the whole yade inside main()
with bunch of loops, if's and switches; it could be done, but it is ugly
and difficult to maintain.

>> (well, let me add 3. Constitutive law should just map strains to
>> stresses. It is not constitutive law's business to compute strains. You
>> can have a look at ef2_Spheres_NormalShear_ElasticFrictionalLaw, this is
>> what I think constitutive law should look like: get strains, compute
>> stresses, apply forces. Finito.)
>>   
> It makes sense, but it means each constitutive law will need its own
> specific geometry engine+dataClass class to compute/store the relevant
> "strain", doesn't it?
> Janek, needs an initial relative orientation, Wenji needed maxUn, JF
> needs the volume of the cell from triangulation (stored at bodies
> level) before he can even define the stiffness, the first Yade
> developers thought everybody would need a contact point defined while
> most of us never used it, etc. And I'm sure the next person
> developping a constitutive law will need another geometrical
> information that nobody is expecting for now. Is it compatible with
> the design you plan?
You can always derive from the geometry class and add your own data
members. But you still have access to inherited data members and
methods, should you need those. maxUn can (should) be stored in
InteractionPhysics, the same for cell volume, no?
>
> Also, when you say "constitutive law should look like: get strains,
> compute stresses, apply forces"(*), I'd say it is perhaps a bit more
> complex.
> Take this equation : Fn = Kn (Un-Un_plastic), obviously an
> elasto-plastic framework.
> What is Un_plastic? Is it an internal parameter of the law or an
> external geometrical information? It is in fact something in between,
> since (1) it  is obviously a kinematic-like variable, but (2) at the
> same time it can't be defined from the current geometry. Un_plastic
> describes an history of the motion _with respect to_ one constitutive
> law. Even with the full record of all previous movements, you can't
> define Un_plastic if you don't know the constitutive law itself.
> Should Un_plastic be in contactPhysics or contactGeometry?
> Do we really need a rule?
Good point. I have SpheresContactGeometry::splitToEpsTMax(Real epsTMax),
which adjusts the original contact points such that |epsT|<epsTMax and
returns increment of eps_plastic. I believe that most constitutive laws
use similar notions of plasticity, so to me it made sense to add this to
SpheresContactGeometry.
> An exemple of this : in creep equations (creep, like plasticity, is
> changing the "neutral" position when it occures), the initial
> orientations (from wich total bending/twisting are defined) are
> updated to reflect creep. The way it is updated is valid only for the
> specific creep law considered here, another constitutive law would
> update the initial configuration differently.
>
> The conclusion is you can't avoid the manipulation/storage of
> geometrical data in some constitutive laws/physicalParameters.
Constitutive law can change geometry (as with the plasticity), in some
special cases. I know that I say something more than just "compute
stresses from strains". I would say that what is useful for multiple
constitutive laws could be defined as SpheresContactGeometry method
(plasticity), otherwise just that single constitutive law can manipulate
that.
> (*) Quick note : I feel strange each time you speak about strain and
> stress in constitutive law, since the specificity of DEM is to
> consider disp./forces instead of strain/stress, usually...
That's just because the models I use are somehow derived from continuum
models. But internally nevertheless, stresses are multiplied by link
cross-sections etc, to it should really be the same after all.
>> What makes you think that contact of 2 spheres is describable just by
>> penetrationDepth and radii?
> Answer 1 : Cundall (1988) ;-)
> Answer 2 : because I need nothing else.
>>  That would be _normal_ contact. For shear,
>> bending and torsion, it gets more complicated, and that's what all the
>> cp1rel, cp2rel (and unrollSpherePtToPlane, for that matter, which
>> computes the shear from the data) are for.
>>   
> The incremental definition of the shear force allows to skip the
> definition and storage of cp1rel/cp2rel/unroll. It works very well and
> it has been used by hundreds of researchers since 1988 (so again, I'd
> like to know how you compared with E.F.).
If you say incremental shear, you still need access to 2 linear
velocities and 2 angular velocities, not just penetrationDepth and radii.

I personally don't trust the incremental code, although I admit that it
can be my personal weakness that I wasn't able to derive it and verify
myself. Besides that, it supposes that the rotation is such small that
displacements can be linearized (if I understand what it does) and I
think that could make for some problems in highly dynamic simulations.

If you think it is realiable and useful, I would be all for extracting
it to a stand-alone function that would get all the parameters it needs
and return what it computes. Or shear strain could be added to
SpheresContactGeometry (or a new class, if you like), but computed by
InteractingSphere2InteractingSphere4SpheresContactGeometry. That would
make more sense to me.
>> The "abusive storage" members were not added by me, except
>> facetContactFace (abusive, because only used for sphere-facet contact).
>> What I agree with is that hasNormalViscosity, NormalViscosity and
>> NormalRelativeViscosity are clearly physical quantitites that should be
>> in interaction physics.   
> So will you move them away at one point?
Oh, I didn't add those (except facetContactFace, which I'd like to keep
until SphresFacetGeometry is done). "svn blame" gives me "richefeu". But
I can remove that, no problem at all ;-)
>> About my future plans: define abstract classes Dem3DofContactGeometry
>> (better name?; for classes normal and shear) and Dem6DofContactGeometry
>> (for classes with normal, shear, bending and torsion strains)
> So, first of all, I would define "DemXDofContactGeometry" with normal,
> Un, and incremental dUt (this code duplicated in different
> constitutive laws), nothing more. This is "normal" DEM, no reason to
> overload with more parameters.
> That is what I am a bit annoyed with the ammount of extra stuff that
> is in SpheresContactGeometry right now.
I put those extra things under condition 'hasShear', they are not
assigned and should "just" take memory if you don't use them (29 bytes,
if I calculate correctly: 5 quaternions, 2 Vector3, 3 Reals). If that is
not OK and if other peoples will not use it, I can move it to a separate
class. I thought it would be a way to clean lot of code elsewhere, though.
>> with
>> virtual methods to get e.g. epsN(), epsT(), regardless of whether it is
>> Sphere-Facet of Sphere-Sphere contact. Mapping sphere-facet to
>> SpheresContactGeometry has not much physical meaning and will not work
>> for more complex things like transferring plastic strain when going over
>> facet boundary to its neighbor.   
> I think I get the idea. I admit I have no clear vision on all this for
> now, I'm just raising a few points.
And I am glad for your comments. Really.
> The possible problem is that people could need to modify the abstract
> classes all the time because they need a geometrical information you
> disregarded
> One simple exemple comes to mind : center of mass of bodies. For
> spheres and simple elastic-frictional contacts, you can compute the
> moment generated by fs as M += radius.n x fs, but for non-sperical
> bodies (snow) you cannot do this. You need to compute M+=(ContPoint -
> CenterMass) x fs.
> In that case, with your approach, IntGeom would need to store not only
> rad1,rad2 but also the centers of mass of both particles.
It is already stored there (iff hasShear): pos1, pos2, ori1, ori2. But I
have yet to see non-spherical particles in yade (if someone is serious
about that, he'd better fix global/local mismatch of inertia in Newton's
law first)
> So, fundamentally, with you first statement in mind : "Constitutive
> law should operate just on interaction geometry and interaction
> physics. If it accesses rootBody->bodies, it is a design problem", all
> information from bodies must be duplicated in interactionGeometry,
> which is a serious design problem too...
I know it is a hassle, it could be that SpheresContactGeometry would
contain shared_ptr<Se3> shared with Body::physicalParameters, but
current serialization wouldn't handle that correctly (it would create 2
copies of the same Se3 instead of the shared one). I still have the
parallel stuff in mind and it is advantageous if the computation is as
local as possible (i.e. no access to rootBody->bodies and such).

>> Then SpheresContactGeometry inherits from Dem3DofContactGeometry, I will
>> create (maybe) also SphereFacetContactGeometry etc, but the constitutive
>> law will operate on Dem3DofContactGeometry and will not care, whether
>> the underlying contact is sphere-sphere, sphere-facet (or sphere-box,
>> whatever).  
> Agreed. The fact that SpheresContactGeometry can be used for others
> (body1,body2) is just a temporary trick.
("temporary" = since as long as I remember...)

Regards, vaclav.