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

Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2021-06-16: #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

 Revision history for this message Liu Changdong (changdong) said on 2021-06-16: #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 on 2021-06-16: #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 on 2021-06-16: #4

Hi Jan:

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 Jan Stránský (honzik) said on 2021-06-16: #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
return predicate(pos,r)
###

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

Cheers
Jan

 Revision history for this message Liu Changdong (changdong) said on 2021-06-16: #6

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

 Revision history for this message Liu Changdong (changdong) said on 2021-06-16: #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
return halfSpace1(pos,r) and halfSpace2(pos,r)

Thanks again!
changdong Liu

 Revision history for this message Liu Changdong (changdong) said on 2021-06-16: #8

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

 Revision history for this message Jan Stránský (honzik) said on 2021-06-16: #9

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

To post a message you must log in.