Yade

How to create a cylinder model with hard and soft interlayer

Asked by Liu Changdong

Hi everyone:

I want to creat a model with hard and soft interlayer like paper [1] Fig 4 . But i do not know how to do that. Can someone help me solve this problem.

Best

changdong Liu

[1] https://www.researchgate.net/publication/338102898_Numerical_simulation_of_uniaxial_compression_tests_on_layered_rock_specimens_using_the_discrete_element_method

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said : 		#1

Hello,

a general approach:
- create "uniform" cylinder
- select "different" (harder or softer) particles
- change their material

> But i do not know how to do that

To get more detailed answer, please read [1] and provide more detailed information what is the actual problem - what you do not know (already creating the cylinder? selecting particles? changing material? ... ?), if you already tried something etc.

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
Liu Changdong (changdong) said : 		#2

Hi Jan：
I'm sorry I didn't make my question clear. My question is as follows:

1. How to select particles representing soft rock (different thickness and dip angle).

2.About changing the material, can i use the following method:

sofeMat=O.materials.append(JCFpmMat(....))
hardMat=O.materials.append(JCFpmMat(....))

pred=pack.inCylinder()
sp=pack.randomDensePack(pred,...)
sp.toSimulation()

'select the soft rock particles'

for i in 'softrock':
   O.bodies[i].mat.cohesion=softMat.cohesion
   ....

Thanks

changdong Liu

Revision history for this message
Jan Stránský (honzik) said : 		#3

> 1. How to select particles representing soft rock (different thickness and dip angle).

e.g.:
###

# return true if z-coordinate is between 1.5 and 2.5
# adjust the condition to your needs
def isSoftRockParticle(b):
    x,y,z = b.state.pos
    return z > 1.5 and z < 2.5

softRockParticles = [b for b in O.bodies if isSoftRockParticle(b)]
###

> 2.About changing the material, can i use the following method:

in general yes, BUT

> for i in 'softrock':
> O.bodies[i].mat.cohesion=softMat.cohesion

does probably something different than you want, better use
###
sp.toSimulation(material=mat1)
for b in softRockParticles:
    b.mat = mat2
###
i.e. change material entirely.

The reason is that material is (usually) shared among (many) particles and changing b.mat.property of one particle applies to the other particles, too (because b.mat is the same Material).
Have a try:
###
mat1 = FrictMat(young=3)
b1 = sphere((0,0,0),1,material=mat1)
b2 = sphere((2,0,0),1,material=mat1)
print(b1.mat.young)
print(b2.mat.young)
b1.mat.young = 7
print(b1.mat.young)
print(b2.mat.young) # !!!
###

compared to

###
mat1 = FrictMat(young=4)
mat2 = FrictMat(young=5)
b1 = sphere((0,0,0),1,material=mat1)
b2 = sphere((2,0,0),1,material=mat1)
print(b1.mat.young)
print(b2.mat.young)
b1.mat = mat2
print(b1.mat.young)
print(b2.mat.young) # OK
###

Cheers
Jan

Revision history for this message
Liu Changdong (changdong) said : 		#4

Hi Jan:

Thank you for your answer.

Problem 2 has been solved. But I have some questions about question 1.

Here are my questions:
According to your answer, I can only create soft rocks with zero dip angle.
How can i create a soft rock layer with 30 dip angle?

Thanks
changdong Liu

Revision history for this message
Best Jan Stránský (honzik) said : 		#5

The posted code was meant to present one option how to do it, with very easy isSoftRockParticle condition.
You can substitute whatever instead of this "z in (1.5;2.5) range".

If your region is bounded by two parallel planes [2], you can easily do something like:
###
normal = Vector3(1,0,4)
def plane(x,y,z,d):
    nx,ny,nz = normal
    return nx*x + ny*y + nz*z + d
def plane1(x,y,z): # equation of "lower" plane
    return plane(x,y,z,2)
def plane2(x,y,z): # equation of "upper" plane
    return plane(x,y,z,3)
def isSoftRockParticle(b):
    x,y,z = b.state.pos
    p1 = plane1(x,y,z)
    p2 = plane2(x,y,z)
    return p1 > 0 and p2 < 0
    # p1 > 0 ... (x,y,z) is "above" the "lower" plane
    # p2 < 0 ... (x,y,z) is "below" the "upper" plane
###

For more complicated situations (or to make your solution easier to extend / maintain / read / ...), you can use Yade predicates [3,4], e.g. pack.inHalfSpace [5]:
###
normal = Vector3(1,0,4)
halfSpace1 = pack.inHalfSpace((0,0,3),+normal) # pack.inHalfSpace(point,normal direction)
halfSpace2 = pack.inHalfSpace((0,0,3),-normal)
predicate = halfSpace1 & halfSpace2 # intersection
def isSoftRockParticle(b):
    pos = b.state.pos
    r = b.shape.radius
    return predicate(pos,r)
###

As a predicate, you can also use more fancy pack.inGtsSurface [6].

Cheers
Jan

[2] https://en.wikipedia.org/wiki/Plane_(geometry)
[3] https://yade-dem.org/doc/user.html#constructive-solid-geometry-csg
[4] https://yade-dem.org/doc/yade.pack.html
[5] https://yade-dem.org/doc/yade.pack.html#yade.pack.inHalfSpace (poor documentation)
[6] https://yade-dem.org/doc/yade.pack.html#yade._packPredicates.inGtsSurface

Revision history for this message
Liu Changdong (changdong) said : 		#6

Hi Jan：
Thanks a lot。
That solves my problem。

Revision history for this message
Liu Changdong (changdong) said : 		#7

Hi Jan：
I have tried the second method you gave me ‘using the pack.inHalfSpace’ . However, the 'predicate(pos,r)' shows an error.

After my trial and error, I found the following method should to be used:

halfSpace1 = pack.inHalfSpace((0.4,0.4,.4),Vector3(1,1,1)) # pack.inHalfSpace(point,normal direction)
halfSpace2 = pack.inHalfSpace((0.6,0.6,.6),Vector3(-1,-1,-1))
#predicate0 = halfSpace1 & halfSpace2 # intersection
def isSoftRockParticle(b):
    pos = b.state.pos
    r = b.shape.radius
    return halfSpace1(pos,r) and halfSpace2(pos,r)

Thanks again!
changdong Liu

Revision history for this message
Liu Changdong (changdong) said : 		#8

Thanks Jan Stránský, that solved my question.

Revision history for this message
Jan Stránský (honzik) said : 		#9

Yes, sorry, I forgot about [7].
Nice you fixed the solution yourself :-)
Cheers
Jan

[7] https://gitlab.com/yade-dev/trunk/-/issues/105