I want to make a cluster of shperes representing a rock

Asked by joao miguel manso on 2012-01-27

Hi,

I want to model a triaxial test with rockfill in it. I want to be able to build clusters of spheres with bonds between them that (representing the rockfill particles), during consolidation and shear phases, could be broken if a threshold value is attained.

I have already tried to write the simplest case I could remember and a the code is as follow:

import math
radius=1

O.bodies.append(utils.geom.facetBox((1,1,1),(5,5,5),wallMask=31))

idConcrete=O.materials.append(CFpmMat(poisson=0.25,young=1e9,density=4800))

O.bodies.append([
 utils.sphere((radius,radius,radius),radius,material=idConcrete),
 utils.sphere((3*radius,radius,radius),radius,material=idConcrete),
        utils.sphere((2*radius,radius+sqrt(3)*radius,radius),radius,material=idConcrete),
 utils.sphere((2*radius,radius+sqrt(3)*radius/3,radius+sqrt(6)*radius*2/3),radius,material=idConcrete),
])

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
 [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom()],
 [Ip2_CFpmMat_CFpmMat_CFpmPhys()],
 [Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM()],
   ),
   GravityEngine(gravity=(0,0,-9.81)),
   NewtonIntegrator(damping=0.4),
]

O.dt=.5e-4*utils.PWaveTimeStep()

When I try to run, it exits yade. I don't know why.
I was able to see the particles falling before, but they are not bonded between them. My aim is to be able to model several rocks (each rock made of bonded spheres) that will break and split if a certain amount of stress is attained.

Thanks,
Joao

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
joao miguel manso
Solved:
2012-01-30
Last query:
2012-01-30
Last reply:
2012-01-30
Christian Jakob (jakob-ifgt) said : #1

I am not sure (because I am also new in Yade comunity), but I think L3Geom is not compatible with ScGeom, so try to use

[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],

instead of

[Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom()],

Good luck!

Jérôme Duriez (jduriez) said : #2

> When I try to run, it exits yade.

Hello, maybe it could be useful that you precise how you did try to run it, and which was the kind of exit ("Normal exit", or "segmentation fault"..) ?

joao miguel manso (jmanso) said : #3

Hi,

In this model, the spheres are supposed to collide with the box and remained bonded.
I already changed IGeomFunctor to [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()] and it actually worked (Yade didn't exit when I tried to run it). However, when the spheres collide with the box, yade exits again with the following message:

terminate called after throwing an instance of 'std::runtime_error'
  what(): Undefined or ambiguous IPhys dispatch for types FrictMat and CFpmMat.
Aborted

Does someone knows what this means?

Jérôme Duriez (jduriez) said : #4

This means, that Yade faced the situation of a contact between the two materials FrictMat and CFpmMat, which it could not handle (to derive the physic of the interaction "IPhys", e.g. normal stifness of the interaction)

Indeed, you took into account the only possibility of contact between CFpmMat and itself with the line : [Ip2_CFpmMat_CFpmMat_CFpmPhys()],

A quick look at your script make me feel that it's because your facets are probably of "FrictMat" (you did not define their mat explicitely..)

joao miguel manso (jmanso) said : #5

jduriez thanks for your answer.
Now that I defined the facets material, it runs, but the spheres are not bonded yet.
How can I bond spheres and define their characteristics in order to simulate a rockfill particle (in this case one made of 4 spheres)?

joao miguel manso (jmanso) said : #6

At this point, my code is like this:

import math
radius=1

idConcrete=O.materials.append(CFpmMat(poisson=0.25,young=1e9,density=4800))

O.bodies.append([
 utils.sphere((radius,radius,radius),radius,material=idConcrete),
 utils.sphere((3*radius,radius,radius),radius,material=idConcrete),
        utils.sphere((2*radius,radius+sqrt(3)*radius,radius),radius,material=idConcrete),
 utils.sphere((2*radius,radius+sqrt(3)*radius/3,radius+sqrt(6)*radius*2/3),radius,material=idConcrete),
])

O.bodies.append(utils.geom.facetBox((1,1,1),(5,5,5),wallMask=31,material=idConcrete))

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
 [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
 [Ip2_CFpmMat_CFpmMat_CFpmPhys()],
 [Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM(label='contactLaw')],
   ),
   GravityEngine(gravity=(0,0,-9.81)),
   NewtonIntegrator(damping=0.4),
]

contactLaw.frictionAngle=.6
contactLaw.FnMax=1e9
contactLaw.FsMax=1e9
contactLaw.isCohesive=True

O.dt=.5e-4*utils.PWaveTimeStep()
O.saveTmp()
#O.step()

I looks fine, but when I check the interactions, in the interactions inspector, none of the contactLaw properties are defined. Does anyone knows why?

Jérôme Duriez (jduriez) said : #7

Because this needs some computations and, hence, at least one iteration to be ran ! ;-)

=> O.run(...) or "play" with the mouse

joao miguel manso (jmanso) said : #8

Ok, what I meant was that, after running the model and the spheres collide, when I check the interactions, in the interactions inspector, none of the contactLaw properties are defined, i.e.: FnMax is 0, as well as FsMax and the isCohesive box is unchecked.
Does anyone knows why?

joao miguel manso (jmanso) said : #9

import math
radius=1

idConcrete=O.materials.append(CFpmMat(poisson=0.25,young=1e9,density=4800))

O.bodies.append([
 utils.sphere((radius,radius,radius),radius,material=idConcrete),
 utils.sphere((3*radius,radius,radius),radius,material=idConcrete),
        utils.sphere((2*radius,radius+sqrt(3)*radius,radius),radius,material=idConcrete),
 utils.sphere((2*radius,radius+sqrt(3)*radius/3,radius+sqrt(6)*radius*2/3),radius,material=idConcrete),
])

O.bodies.append(utils.geom.facetBox((1,1,1),(5,5,5),wallMask=31,material=idConcrete))

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
 [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
 [Ip2_CFpmMat_CFpmMat_CFpmPhys(cohesion=1000,tensileStrength=1e9,Alpha=1)],
 [Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM(label='contactLaw')],
   ),
   GravityEngine(gravity=(0,0,-9.81)),
   NewtonIntegrator(damping=0.4),
]

O.dt=.5e-4*utils.PWaveTimeStep()
O.saveTmp()