Writing a new constitutive law

Hi,

why can it be painful? Did you try it? Yade is actually designed in a way that new contact laws can easily be added and it is open source so you are in control of everything. The only thing you might need is a little bit of C++ and Python knowledge.

You should first read the Programmer's manual [1] and then have a look at the source code in pkg/dem/ where all the files defining the contact laws are. Choose the contact law which is most similar to your new one (if you tell us more about your new contact law we might be able to give you some more directions) and start from there, i.e. copy and rename the relevant files and start to modify them. Make sure you use the right namings for your classes as indicated in the Programmer's manual.

Good luck,
Klaus

[1] https://yade-dem.org/doc/prog.html

Hello,
If you have the equations written down on a paper, you can show them here so we can anticipate the possible difficulties and advise on reusing code fragments or functors.
B

Thanks guys for the answers. I start reading the programmer manual and I will get back to you if I get stuck somewhere.
I really appreciate your helps.

cheers,
Behzad

Hey,

I know python and I'm learning C++. I do not know the steps I need to do to include my new things (can be a contact law or something else) in Yade source code.

There're some issues in Yade for me now;

1-We assign material properties to "Box"es and so they've got weight and they follow Newton's second law. But, momentum is not considered in their movement law!

2-We assign material properties to Boxes but why not to facets?

The idea behind all these is to import some aggregates which we usually define as clumps in the form of facets and assign material properties to them (only for the movement, they will be considered as rigid bodies). In models with thousands of clumps this will reduce the number of elements and there will be a computational gain. The other phases of the material will be defined with spheres as usual.

So, what I'd like to do is to be able to import some shapes as gts and assign a density to them and include it in Newton integration engine for movement calculations.
There might be a need for new contact detection definitions like facet_to_box.

Hi Behzad,

> I know python and I'm learning C++. I do not know the steps I need to do
> to include my new things (can be a contact law or something else) in
> Yade source code.
>

For the beginning, define a small task for yourself and share it with us,
we can guide you. The better and more specific is the description you give
us the better can be our understanding of your problem and our help :-)

>
> There're some issues in Yade for me now;
>
> 1-We assign material properties to "Box"es and so they've got weight and
> they follow Newton's second law. But, momentum is not considered in
> their movement law!
>

I am not sure what you mean by this statement.. The second Newton's law is
momentum balance, so the momentum of bodies is essentially present in the
solution.

>
> 2-We assign material properties to Boxes but why not to facets?
>
>
material properties (instance of Material class) is assigned to any body.
Facets are non-dynamic by default, do you mean this? It is because people
usually use them as boundary with prescribed position. But for interactions
evaluation, they do have material properties.

> The idea behind all these is to import some aggregates which we usually
> define as clumps in the form of facets and assign material properties to
> them (only for the movement, they will be considered as rigid bodies).
> In models with thousands of clumps this will reduce the number of
> elements and there will be a computational gain. The other phases of the
> material will be defined with spheres as usual.
>

do you mean something like polyhedrons?

>
> So, what I'd like to do is to be able to import some shapes as gts and
> assign a density to them and include it in Newton integration engine for
> movement calculations.
> There might be a need for new contact detection definitions like
> facet_to_box.
>

One way is to use existing implementation of polyhedrons for this task (can
be used for both "gts shapes" and boxes). Depending on what you will write
us about what you want and expect we can suggest further steps.

Sorry if I misunderstood something in your question, it's quite late ;-)

cheers
Jan

Hi

About your question "1-", I did not check what you said, but I would not be affraid of such things.

Until now, boxes (as facets) were, I guess, always used as boundaries elements whose velocities are directly prescribed. Their movements thus do not depend on sum of acting forces. (By the way, for this reason, I think that in fact 2d Newton's law is not applied at all to boxes).
Have in mind that NewtonIntegrator changes the position of any body according to its velocity. Concerning velocity, depending on the dynamic feature / the blocked DOF, the velocity is either computed according to acceleration = sum of forces / mass, or constant (equal to 0 or the user-defined value )

About 2-, I do not use facets, but I also think that they also have material properties. "Facet" is just a possible "shape" for "Body", that has always a "Material"

as for Jan;

What I mean about Boxes is that, we have the option of "fixed=False" which allows it to move and since it has already density assigned, it can and it does move because of gravity acceleration. However, I think it has only 3 DOF's.
I wanted to see if we can change facets to dynamic bodies.
Anyway, this is not important if as you guys said I can use Polyhedra for the purpose of grains.

Yes, what I need to to is import gts or stl into Yade and treat them the same way as polyhedra. The problem with current implementation of polyhedra is that we cannot have that much control on its shape. But I have already certain irregular-shape grains as stl.

as for jduriez;

The reason I wanted to assign properties such as density to facets is that (as I mentioned above for Jan) to import gts as facets and treat them as polyhedra of clumps! So I wanted to make them Dynamic bodies.

I think, according to all I can conclude that the best way to deal with the issue is using polyhedra. To do so I will need;

1- find a way to import gts/stl and define it as polyhedra in Yade
2- Define a polyhedra-sphere contact law

Jan;

To start writing and compiling a new constitutive law, perhaps it's not a bad idea if I first have a look on the simple already implemented laws.
for example where I can find the C++ files of a simple law of frictMat=frictMat contact?
Let's say I modify this, then how do I compile it?

Thanks a lot,

Behzad

Hi Behzad,

> What I mean about Boxes is that, we have the option of "fixed=False" which
> allows it to move and since it has already density assigned,

to be precise, the Body (its State) has mass (not density) asign :-) which
is usually computed according to Material density and Shape volume, but you
can set any value of mass if you want.

> it can and it does move because of gravity acceleration. However, I think
> it has only 3 DOF's.
>

Apart from mass, State also has intertia parameter, playing role in
evaluation of rotational DOFs. All shapes in Yade has AFAIK 6 DOFs

> I wanted to see if we can change facets to dynamic bodies.
> Anyway, this is not important if as you guys said I can use Polyhedra for
> the purpose of grains.
>

Facets as planar elements has zero volume and therefore would have zero
mass and inertia if computed from density. However, you can set it manually
to make it work:

f = facet(...)
f.dynamic = True
f.state.mass = someMass
f.state.inertia = someInertia
f.state.ori = someOri # as inertia is Vector3 of principal values, you have
to set initial orientation to make is consistent
O.bodies.append(f)

now the facet would behave dynamically

>
> Yes, what I need to to is import gts or stl into Yade and treat them the
> same way as polyhedra. The problem with current implementation of
> polyhedra is that we cannot have that much control on its shape. But I
> have already certain irregular-shape grains as stl.
>

you can import both gts and stl files to create facets (see corresponding
files yade/lib/import/STLReader.hpp and yade/py/ymport.py). It should not
be very difficult to implement it also for polyhedrons. In the case you
already know vertices of your polyhedron, you can have full control over
its shape. Have a look as spheres and facets are constructed in utils
module [1,2]. In your case it would be analogous to facet, but with more
vertices The only limitation is that in current implementation the
polyhedron has to be convex.

>
> as for jduriez;
>
> The reason I wanted to assign properties such as density to facets is
> that (as I mentioned above for Jan) to import gts as facets and treat
> them as polyhedra of clumps! So I wanted to make them Dynamic bodies.
>

Anyway I think it should be possible to construct your clump from facets
and then assign proper mass and inertia to it. Then it could work for
interactions with speheres, but there is no facet-facet interaction
implemented yet :-) So your facet clumps could not interact with each
other.. This paragraph is only my opinion as I have never tried this. Also
in clump source file there are several comments that clump members should
be spheres, but I did not study them in detail...

>
>
> I think, according to all I can conclude that the best way to deal with
> the issue is using polyhedra. To do so I will need;
>
> 1- find a way to import gts/stl and define it as polyhedra in Yade
> 2- Define a polyhedra-sphere contact law
>

You have two options:
1) use polyhedrons, but only polyhedron-polyhedron interaction is
implemented. For polyhedron-anyOtherShape you would have to implement a new
features
2) use facets (well, firstly ask another question if it is
implemented/possible to implement clump made of facets) and implement new
facet-facet interactions.

In any way, importing gts/stl file is not the difficult part of this task
:-)

cheers
Jan

[1]
http://bazaar.launchpad.net/~yade-pkg/yade/git-trunk/view/head:/py/utils.py#L142
[2]
http://bazaar.launchpad.net/~yade-pkg/yade/git-trunk/view/head:/py/utils.py#L120

1 Compile (https://yade-dem.org/doc/installation.html)
2. mess up a contact law
3. compile again (make install)

If you would like to make this new law public, which would be good for everyone but especially for you, come back to us for comments and advices.

Check the kedevelop page at yade's wiki: https://yade-dem.org/wiki/Using_kdevelop4
It gives plenty ways to find were a law is defined (ctrl+alt+f, "quick open" feature, ...)

Hi Behzad,

> Jan;
>

this can be any developer :-)

> To start writing and compiling a new constitutive law, perhaps it's not a
> bad idea if I first have a look on the simple already implemented laws.
>

yes :-)

> for example where I can find the C++ files of a simple law of
> frictMat=frictMat contact?
>

the answer to your question is not so easy.. have a look at [1,2] chapters
in user's manual. For the standard basic case:
- when two bodies newly interct:
- Firstly Ig2 (e.g. Ig2_Sphere_Sphere_ScGeom) functor is called creating
geometric properties of the interaction from bodies' Shapes.
- these properties are stored in IGeom instace, belonging to interaction
(in this case ScGeom)
- then Ip2 functor is called, in this case Ip2_FrictMat_FrictMat_FrictPhys.
It take the two bodies' Materials and create new IPhys (specifically
FrictPhys) instance.
- IPhys stores variables specific for the interaction physics (could be
e.g. material parameter like internal friction angle - same for all
interactions - but also internal variables like plastic strain - different
for every interaction -etc.)
- for forces evaluation, you need Law2 functor (in this case
Law2_ScGeom_FrictPhys_CundallStrack) combining IGeom and IPhys of the
interaction and computing forces, that are then used by NewtonIntegrator

Files related to the example:
pkg/common/Sphere.xpp
pkg/dem/Ig2_Sphere_Sphere_ScGeom.xpp
pkg/dem/ScGeom.xpp
pkg/common/ElastMat.xpp
pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.xpp
pkg/dem/FrictPhys.xpp
pkg/dem/ElasticContactLaw.xpp

Let's say I modify this, then how do I compile it?
>

See [3]. After a modification, recompile yade by 'make' and 'make install'
commands

in case of any problem, write us, we will try to help
cheers
Jan

[1] https://yade-dem.org/doc/user.html#creating-interactions
[2] https://yade-dem.org/doc/user.html#base-engines
[3] https://yade-dem.org/doc/installation.html#source-code

Hi

I checked out the files like;

Files related to the example:
pkg/common/Sphere.xpp
pkg/dem/Ig2_Sphere_Sphere_ScGeom.xpp
pkg/dem/ScGeom.xpp
pkg/common/ElastMat.xpp
pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.xpp
pkg/dem/FrictPhys.xpp
pkg/dem/ElasticContactLaw.xpp

Alright, now I know how can I compile and use the modifications I might do.

But, for example in which of those mentioned files are the simple equation of f=k*deltaU?

Imagine I want to double the contact force in the simple elastic-elastic contact! It means let say instead of f=k*deltaU I wanna put f=2*k*deltaU.

Hi Behzad,

> But, for example in which of those mentioned files are the simple
> equation of f=k*deltaU?
>
> Imagine I want to double the contact force in the simple elastic-elastic
> contact! It means let say instead of f=k*deltaU I wanna put
> f=2*k*deltaU.

what you want is evaluation of forces, i.e. what you need is Law2 (see also
previous email or user's manual). Specifically
pkg/dem/ElasticContactLaw.cpp, function
void Law2_ScGeom_FrictPhys_CundallStrack::go

cheers
Jan

Alright, from what I understood so far;

Let's say I want to create a new constitutive law in which the contacts between the spheres are governed by a Kelvin model(dashpot in parallel with a spring) in series with another spring.

I this case I need a new InteractionPhysics.
The contact law will have three parameters two spring stiffnesses and a dashpot viscosity. It is supposed to be a bulk material ans so will need bonds between the spheres.

Do I need to define a new class of material?
Or we define a functor which creates the required InteractionPhysics?
Or both, even?!

We have already CpFrictMat which holds normal and shear cohesion. Can I use that?

@Jan: "Anyway I think it should be possible to construct your clump from facets
and then assign proper mass and inertia to it" From some studies of sphere drop simulations I'd like to note that contact force calculation at "edges" (interactin of facets) is not correct, at least in the case of a contacting sphere [1]. We had a discussion on that some time ago. In the current facet implementation edges are not defined. That is why at the intersection of multiple facets sitffness is too high (sum of the stiffnesses of all intersecting facets if I recall correctly).
As a consequence there will be an increasing error in contact forces with increasing number of facet-intersections. Imagine a sphere approximated by facets. If you use lets say 48 instead of 12 facets for approximation of the sphere you will have way more intersections of facets (edges).

behzad: "It is supposed to be a bulk material ans so will need bonds between the spheres."
why do you need bonds between spheres for bulk material? What kind of material is it? If you consider a packing of spheres that is deposited by gravity there don't need to be bonds.

[1] https://answers.launchpad.net/yade/+question/229901

The material is not granular and it's pretty much like asphalt binder. So, we need bonds between the spheres.

Hi,

the new constitutive law I need to implement is that of Burger's model. This model is composed of a Kelvin model (dashpot in parallel of a spring) in series of a Maxwell model (dashpot in series with a spring). Thus, it has four parameters; two viscosity of the dashpots and two stiffness for two springs.

The deformation equation of this model is given by;

http://books.google.ca/books?id=E1ZkAkDjbuIC&pg=PA183&dq=strain+function+burgers+model&hl=en&sa=X&ei=wsl0U4uvFcyBqgbbl4H4Bw&ved=0CDQQ6AEwAQ#v=onepage&q=strain%20function%20burgers%20model&f=false

I don't know how should I start writing this in C++ for yade!

cheers,
Behzad

Hi Behzad,

I think that you will need new Material, IPhys, Ip2 and Law classes. Have a
look at FrictPhys related classes for an inspiration (see one of previous
answers). Try to understand the relations between those classes. If
everything is clear to you, the extansion to Burger's model should not be a
big deal :-)
cheers
Jan

2014-05-15 16:11 GMT+02:00 behzad <email address hidden>:

> Question #247006 on Yade changed:
> https://answers.launchpad.net/yade/+question/247006
>
> Status: Answered => Open
>
> behzad is still having a problem:
> Hi,
>
> the new constitutive law I need to implement is that of Burger's model.
> This model is composed of a Kelvin model (dashpot in parallel of a
> spring) in series of a Maxwell model (dashpot in series with a spring).
> Thus, it has four parameters; two viscosity of the dashpots and two
> stiffness for two springs.
>
> The deformation equation of this model is given by;
>
>
> http://books.google.ca/books?id=E1ZkAkDjbuIC&pg=PA183&dq=strain+function+burgers+model&hl=en&sa=X&ei=wsl0U4uvFcyBqgbbl4H4Bw&ved=0CDQQ6AEwAQ#v=onepage&q=strain%20function%20burgers%20model&f=false
>
> I don't know how should I start writing this in C++ for yade!
>
> cheers,
> Behzad
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

There is an almost-Burger model implemented in Law2_ScGeom6D_CohFrictPhys_CohesionMoment, actually a "Standard Linear Solid model" [1]. It could help for a starting point (grep "creep" in Yade's code).
I don't know if you are interested in "creeping" the same components though (currently the shear force and twisting moment). But if you want the burger model on the normal force it is very much easier.

I know people in Knoxville already implemented a burger model (I've seen a PhD thesis about this somewhere, it is in Yade's publications IIRC). It is a pitty that you need to start again.

[1] http://en.wikipedia.org/wiki/Standard_linear_solid_model

Hi Bruno,

Yeah, that would be a good point to start if I use (and extend to the Burger's) "Standard Linear Solid model".

which file is describing this contact law?

in my trunk/pkg/dem I have only :

Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp

Hi,

Law2_ScGeom6D_CohFrictPhys_CohesionMoment is in pkg/dem/CohesiveFrictionalContactLaw.* (hpp or cpp)
Searching the string Law2_ScGeom6D_CohFrictPhys_CohesionMoment, through kdevelop for instance, helps you to locate it.

Jérôme

Bruno,

what you mentioned, Standard linear Solid, contains two springs and a dashpot.
How the parameters of this model (dashpot viscosities and spring stiffness) are defined?

What I see in yade is;

list(yade.wrapper.Law2_ScGeom6D_CohFrictPhys_CohesionMoment.__dict__.keys())

 -> [4]:
['shearElastEnergy',
 '__module__',
 'shear_creep',
 '__instance_size__',
 '__reduce__',
 'twist_creep',
 'creep_viscosity',
 'always_use_moment_law',
 'useIncrementalForm',
 'normElastEnergy',
 'neverErase',
 '__doc__',
 '__init__']

and other questions;

1- This elements act both in shear and normal direction? How we specify this?
2- What about the cohesion in normal and shear direction? How we define it?

thanks

Hi,

Note that you can also give a look (did you ?) to the doc of this law in [1] to have the list of the parameters.

Normal and shear cohesions are stored in CohFrictPhys (the physics of the interaction, as usual), in [2] and [3], and they depend on CohFrictMat (material) parameters. You have some explanation in the doc of the corresponding Ip2 [4], which does the computations of physical "parameters" of the interaction . If you want even more details about how different "parameters" are computed, you have probably to look in the source code.

This design is always the same in Yade.

Jérôme

[1] : https://www.yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Law2_ScGeom6D_CohFrictPhys_CohesionMoment
[2] : https://www.yade-dem.org/doc/yade.wrapper.html#yade.wrapper.CohFrictPhys.normalAdhesion
[3] : https://www.yade-dem.org/doc/yade.wrapper.html#yade.wrapper.CohFrictPhys.shearAdhesion
[4] : https://www.yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Ip2_CohFrictMat_CohFrictMat_CohFrictPhys

Mmmmh... I realize that the version in trunk is not the last one I used.
Presently it is only a Maxwell model (actually it is the most difficult
part). I'll see if I can find the last version somewhere.
Bruno

Bruno,

Are you talking about "Law2_ScGeom6D_CohFrictPhys_CohesionMoment"?
You mean this is currently a Maxwell model istead of Standard Leaner Solid?

Yes, the github version has a Maxwell model on shear force and twist moment.
B

I was confused, sorry. My previous statements are correct about
Law2_ScGeom6D_CohFrictPhys_CohesionMoment.
But there actually is also a model with two springs and one damper
(Standard Linear Solid [1]) on the shear force. It is in
Law2_ScGeom_ViscoFrictPhys_CundallStrack, and "creepStiffness" is the E2
in [1].
This law has been used by a master student but it is still a bit
experimental. You can improve it if you like.

Bruno

[1] http://en.wikipedia.org/wiki/File:SLS.svg

I'm sorry, but two questions to avoid misunderstandings;

1- Law2_ScGeom6D_CohFrictPhys_CohesionMoment is a Maxwell model and Law2_ScGeom_ViscoFrictPhys_CundallStrack is the Standard Linear Solid, right?

2- This question concerns both models; We will have the same elements in normal and shear directions, right? I mean in Standard Linear model we have 2 springs and a dashpot in normal direction, and 2 springs and one dashpot in shear direction.
Is that correct? It's the case for the models in PFC actually. and the ratio of normal to shear stiffnesses comes from the Poisson ratio.

Thanks

1. yes
2. no. Only shear forces and twisting moment would have a viscous
behaviour at the moment.

Thank you!