Ask help for the materials setting

Asked by yang yi on 2019-12-25

I hope to simulate the coal seam collapsing, but my major is computer and control. And I tried several different material about the particle. such as FricMat, CpmMat, ViscElMat.... But all of them has the following problme
(1) The particle will bounce if it collide the wall of boundary.
(2) My friend to me to increase the damping of NewtonIntegrator, but it desnot work. and if I increase the greavity, the particle will through the wall.

my engines is :
         ForceResetter()
         InsertionSortCollider()

So, please kindly tell me, which kind of material I should to use. and how to set the parameters that the particle will not bounce.
Thank you very much

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
yang yi
Solved:
2020-01-07
Last query:
2020-01-07
Last reply:
yang yi (yangyizaixian) said : #1

Sorry, maybe my description is not clear.
My aim is to establish a environment to simulation the process of coal and rock falling. The upper layer is rock, then a coal layer. The particles falls by gravity. My question is

(1) The particles will damp, but I want that the coal and rock falls down and stays on the ground, or the bounce is very small. So I do not know which kind of material I should to use. I tried the FricMat, CpmMat, ViscElMa...
      Today, I try the "CohesiveDeformableElementMaterial", but there is a error that
      "Undefined or ambiguous IPhys dispatch for types CohesiveDeformableElementMaterial and CohesiveDeformableElementMaterial." So, here is subproblem: If the "CohesiveDeformableElementMaterial" is right, which kind of Engine should I used??

(2) Sometimes , the particles will through the wall.

Please kindly help me.

My code
#!/usr/bin/env python
#  encoding: utf-8
from __future__ import print_function
import sys
sys.path.append('~/PycharmProjects/Yade20191229/')
# from yadeImport import *
import math
import numpy as np

o = Omega()
tc = 0.01
o.dt = 0.02 * tc
numWinds = 5
VelocityWind = 20
WinAction = np.zeros(5, dtype=int)
## =========== Environment ============##
widthSpace = 10.5 ## width is the distance from the front coal wall to the back
widthHydr = 1.5
lengthSpace = widthHydr * numWinds ## length is the width of 5 windows
highCoal = 2

highRock = 1
highHydr = 3.8
highUnderground = 12
highSpace = highUnderground + highHydr + highCoal + highRock
radiusCoal = 0.15
radiusRock = 0.2
CheckThick = 1

colorCoal = (0, 1, 0)
colorRock = (1, 0, 0)
colorNote = (0, 0, 1)
colorHydr = (1, 0, 1)
colorShield = (1, 0, 1)
colorWind = (0, 1, 1)
colorGround = (1, 1, 0)
angleShield = 50 * 3.1415926 / 180
angleSwingPositive = 15 * 3.1415926 / 180
angleSwingNegtive = 40 * 3.1415926 / 180
lengthShield = 3
lengthTail = 2
windUpperBoundary = 0.9
windLowerBoundary = 0.5019
stateUpperBoundary = highHydr
stateLowerBoundary = windLowerBoundary
positionShield = []
positionWind = []
positionTopBeam = []

shield_y_0 = lengthTail * math.cos(angleShield - angleSwingPositive)
shield_y_1 = shield_y_0 + lengthShield * math.cos(angleShield)
shield_z_0 = highHydr - lengthShield * math.sin(angleShield)
shield_z_1 = highHydr #
wind_y_0 = lengthTail * (math.cos(angleShield - angleSwingPositive) - math.cos(angleShield))
wind_y_1 = shield_y_0
wind_z_0 = highHydr - (lengthShield + lengthTail) * math.sin(angleShield)
wind_z_1 = highHydr - lengthShield * math.sin(angleShield)
topBeam_y_0 = lengthTail * math.cos(angleShield - angleSwingPositive) + lengthShield * math.cos(angleShield)
topBean_y_1 = widthSpace

##=============================###
### problem 1 #####
matCoal = O.materials.append(ViscElMat(frictionAngle=0.5, mR=0.05, mRtype=1, density=20000, tc=0.001, en=0.7, et=0.7))
matRock = O.materials.append(ViscElMat(frictionAngle=0.5, mR=0.05, mRtype=1, density=20000, tc=0.001, en=0.7, et=0.7))
mat1 = O.materials.append(FrictMat(frictionAngle=0.1, density=2000, young=2e10))
def Boundary():
    boundary = O.bodies.append(geom.facetBox((lengthSpace / 2, widthSpace / 2, highSpace / 2 - highUnderground),
                                             (lengthSpace / 2, widthSpace / 2, highSpace / 2), wallMask=63))
def Ground():
    positionGround = [
        Vector3(0,
                lengthTail * (math.cos(angleShield - angleSwingPositive) - math.cos(angleShield + angleSwingNegtive)),
                0),
        Vector3(lengthSpace,
                lengthTail * (math.cos(angleShield - angleSwingPositive) - math.cos(angleShield + angleSwingNegtive)),
                0),
        Vector3(lengthSpace, widthSpace, 0),
        Vector3(0, widthSpace, 0)
    ]
    Ground = pack.sweptPolylines2gtsSurface([positionGround], capStart=True, capEnd=True)
    IDGround = O.bodies.append(pack.gtsSurface2Facets(Ground, color=colorGround))

def HydraulicSupport():
    for i in range(0, numWinds):
        temp = [Vector3(widthHydr * i, shield_y_0, shield_z_0),
                Vector3(widthHydr * (i + 1), shield_y_0, shield_z_0),
                Vector3(widthHydr * (i + 1), shield_y_1, shield_z_1),
                Vector3(widthHydr * i, shield_y_1, shield_z_1)
                ]
        positionShield.append(temp)
        temp = [
            Vector3(widthHydr * i, wind_y_0, wind_z_0),
            Vector3(widthHydr * (i + 1), wind_y_0, wind_z_0),
            Vector3(widthHydr * (i + 1), wind_y_1, wind_z_1),
            Vector3(widthHydr * i, wind_y_1, wind_z_1)
        ]
        positionWind.append(temp)
        temp = [
            Vector3(widthHydr * i, topBeam_y_0, highHydr),
            Vector3(widthHydr * (i + 1), topBeam_y_0, highHydr),
            Vector3(widthHydr * (i + 1), topBean_y_1, highHydr),
            Vector3(widthHydr * i, topBean_y_1, highHydr)
        ]
        positionTopBeam.append(temp)
    # kwBoxes={'color':[1,0,0],'wire':False,'dynamic':False,'material':0}
    # vibrationRotationPlate = O.bodies.append(geom.facetBox((-15,5,-5),(2,2,2),wallMask=16,**kwBoxes))
    Shield1 = pack.sweptPolylines2gtsSurface([positionShield[0]], capStart=True, capEnd=True)
    Shield2 = pack.sweptPolylines2gtsSurface([positionShield[1]], capStart=True, capEnd=True)
    Shield3 = pack.sweptPolylines2gtsSurface([positionShield[2]], capStart=True, capEnd=True)
    Shield4 = pack.sweptPolylines2gtsSurface([positionShield[3]], capStart=True, capEnd=True)
    Shield5 = pack.sweptPolylines2gtsSurface([positionShield[4]], capStart=True, capEnd=True)

    IDShield1 = O.bodies.append(pack.gtsSurface2Facets(Shield1, color=colorShield))
    IDShield2 = O.bodies.append(pack.gtsSurface2Facets(Shield2, color=colorShield))
    IDShield3 = O.bodies.append(pack.gtsSurface2Facets(Shield3, color=colorShield))
    IDShield4 = O.bodies.append(pack.gtsSurface2Facets(Shield4, color=colorShield))
    IDShield5 = O.bodies.append(pack.gtsSurface2Facets(Shield5, color=colorShield))

    Wind1 = pack.sweptPolylines2gtsSurface([positionWind[0]], capStart=True, capEnd=True)
    Wind2 = pack.sweptPolylines2gtsSurface([positionWind[1]], capStart=True, capEnd=True)
    Wind3 = pack.sweptPolylines2gtsSurface([positionWind[2]], capStart=True, capEnd=True)
    Wind4 = pack.sweptPolylines2gtsSurface([positionWind[3]], capStart=True, capEnd=True)
    Wind5 = pack.sweptPolylines2gtsSurface([positionWind[4]], capStart=True, capEnd=True)

    IDWind1 = O.bodies.append(pack.gtsSurface2Facets(Wind1, color=colorWind))
    IDWind2 = O.bodies.append(pack.gtsSurface2Facets(Wind2, color=colorWind))
    IDWind3 = O.bodies.append(pack.gtsSurface2Facets(Wind3, color=colorWind))
    IDWind4 = O.bodies.append(pack.gtsSurface2Facets(Wind4, color=colorWind))
    IDWind5 = O.bodies.append(pack.gtsSurface2Facets(Wind5, color=colorWind))

    TopBeam1 = pack.sweptPolylines2gtsSurface([positionTopBeam[0]], capStart=True, capEnd=True)
    TopBeam2 = pack.sweptPolylines2gtsSurface([positionTopBeam[1]], capStart=True, capEnd=True)
    TopBeam3 = pack.sweptPolylines2gtsSurface([positionTopBeam[2]], capStart=True, capEnd=True)
    TopBeam4 = pack.sweptPolylines2gtsSurface([positionTopBeam[3]], capStart=True, capEnd=True)
    TopBeam5 = pack.sweptPolylines2gtsSurface([positionTopBeam[4]], capStart=True, capEnd=True)

    IDTopBeam1 = O.bodies.append(pack.gtsSurface2Facets(TopBeam1, color=colorShield))
    IDTopBeam2 = O.bodies.append(pack.gtsSurface2Facets(TopBeam2, color=colorShield))
    IDTopBeam3 = O.bodies.append(pack.gtsSurface2Facets(TopBeam3, color=colorShield))
    IDTopBeam4 = O.bodies.append(pack.gtsSurface2Facets(TopBeam4, color=colorShield))
    IDTopBeam5 = O.bodies.append(pack.gtsSurface2Facets(TopBeam5, color=colorShield))
    IDWind = [IDWind1, IDWind2, IDWind3, IDWind4, IDWind5]
    return IDWind

def CoalLayer():
    ## establish coal layer
    IDCoal = O.bodies.append(
        pack.regularHexa(inAlignedBox((0, 0, highHydr), (lengthSpace, widthSpace, highCoal + highHydr)),
                         radius=radiusCoal, gap=0, color=colorCoal, material=matCoal))
    return IDCoal

def RockLayer():
    ## estatblish the rock layer
    IDRock = O.bodies.append(
        pack.regularHexa(
            inAlignedBox((0, 0, highHydr + highCoal - radiusCoal * 2),
                         (lengthSpace, widthSpace,
                          highCoal + highHydr + highRock)),
            radius=radiusRock, gap=0, color=colorRock, material=matRock))
    return IDRock

def WindowsAction(IDWind):
    global WinAction
    RotationW = [RotationW1, RotationW2, RotationW3, RotationW4, RotationW5]
    for nW in range(0, numWinds):
        ## action
        Pos_z = sum(O.bodies[i].state.pos[2] for i in IDWind[nW]) / len(IDWind[nW])
        if WinAction[nW] == 0:
            RotationW[nW].angularVelocity = -VelocityWind
        else:
            RotationW[nW].angularVelocity = VelocityWind

        NegtiveStop = (RotationW[nW].angularVelocity > 0) & (Pos_z <= windLowerBoundary)
        PostiveStop = (RotationW[nW].angularVelocity < 0) & (Pos_z >= windUpperBoundary)

        if NegtiveStop | PostiveStop:
            RotationW[nW].angularVelocity = 0

Ground()
Boundary()
IDWind = HydraulicSupport()
IDCoal = CoalLayer()
IDRock = RockLayer()
nCoal = IDCoal
nRock = IDRock
O.engines = [
    ForceResetter(),
    ##======================================##
    ## problem 2 ######
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]),
    NewtonIntegrator(gravity=(0, 0, -200), damping=0.001, label='down'),

    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[0], zeroPoint=positionWind[0][2],
                   label='RotationW1'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[1], zeroPoint=positionWind[1][2],
                   label='RotationW2'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[2], zeroPoint=positionWind[2][2],
                   label='RotationW3'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[3], zeroPoint=positionWind[3][2],
                   label='RotationW4'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[4], zeroPoint=positionWind[4][2],
                   label='RotationW5'),
    PyRunner(command="WindowsAction(IDWind)", iterPeriod=5),

]
# o.run(20000000,True)

Jan Stránský (honzik) said : #2

Hello,

search for "coefficient of restitution" or something like that in question or manual.

another option is to develop / implement a new model suitable for your needs :-) which would take significant amount of time, however..

cheers
Jan

yang yi (yangyizaixian) said : #3

Hi, Jan Stránský (honzik), thank you for your responce. I have find the answer from the question.

Jan Stránský (honzik) said : #4

> I have find the answer from the question.

for future reference (if somebody else comes here with same/similar problem), could you please describe the answer?
thanks
Jan