yade-users team mailing list archive

[Fu zuoguang <question258763@xxxxxxxxxxxxxxxxxxxxx>]

Question #258763 on Yade changed:
https://answers.launchpad.net/yade/+question/258763

Fu zuoguang gave more information on the question:
# unicode: UTF-8
from yade import pack,os,utils
import math

psd_type = input('Please input the type of psd data, 0 or 1: ')
contact_type = input('Shall we consider the relative particle rolling, 0 or 1: ')

# Define two type of psd data for testing
psd_0_data = [[0.3600,0.4000,0.4400,0.4800,0.5200,0.5600,0.6000,0.6400,0.6800,0.7200,0.7600,0.8000,0.8400,0.8800,0.9200,0.9600,1.0000],              [0.0000,0.0004,0.0027,0.0114,0.0384,0.1032,0.2242,0.4002,0.5997,0.7757,0.8967,0.9615,0.9885,0.9972,0.9995,0.9999,1.0000]]

psd_1_data =
[[0.8400,0.8500,0.8600,0.8700,0.8800,0.8900,0.9000,0.9100,0.9200,0.9300,0.9400,0.9500,0.9600,0.9700,0.9800,0.9900,1.0000],
[0.0000,0.0004,0.0027,0.0114,0.0384,0.1032,0.2242,0.4002,0.5997,0.7757,0.8967,0.9615,0.9885,0.9972,0.9995,0.9999,1.0000]]

# Define two type of particle materials, one is for common frictmat and the other is for cohfrictmat. 
particles_mat_0 = O.materials.append(FrictMat( young = 1e9,poisson = 0.5,frictionAngle = 0.59,density  = 2650))

particles_mat_1 = O.materials.append(CohFrictMat( young = 1e9,poisson = 0.5,frictionAngle = 0.59,density = 2650,
                                                  isCohesive = False,momentRotationLaw = True,alphaKr = 1.0,etaRoll = 1.0))
# Define the walls
walls_mat = O.materials.append(FrictMat(young = 1e9,poisson = 0.5,frictionAngle = 0,density=0))
mincorner = (0, 0, 0)
maxcorner = (0.2, 0.6, 0.2)
walls_width = 1e-9
genmincorner = [mincorner[0]+walls_width,mincorner[1]+walls_width,mincorner[2]+walls_width]
genmaxcorner = [maxcorner[0]-walls_width,maxcorner[1]-walls_width,maxcorner[2]-walls_width]
walls = aabbWalls([mincorner,maxcorner],thickness=walls_width,material=walls_mat)
wallids = O.bodies.append(walls)

# Define the boundary condition for running 
consolidation = TriaxialStressController(
    wall_left_id=wallids[0],wall_right_id=wallids[1],
    wall_bottom_id=wallids[2],wall_top_id=wallids[3],
    wall_back_id=wallids[4],wall_front_id=wallids[5],
    wall_front_activated = False,wall_back_activated = False,
    maxMultiplier = 1.01,           
    finalMaxMultiplier = 1.001,     
    thickness = walls_width,
    internalCompaction = True,        
    stressMask = 7,
        goal1 = - 100000,  
        goal2 = - 100000,  
        goal3 = - 100000,  
)

# particle generation with a given psd
sp0 = pack.SpherePack();
if psd_type == 0:
    sp0.makeCloud(genmincorner,genmaxcorner,num = 3000,
                  psdSizes = psd_0_data[0],psdCumm = psd_0_data[1],
                  distributeMass = False)
elif psd_type == 1:
    sp0.makeCloud(genmincorner,genmaxcorner,num = 3000,
                  psdSizes = psd_1_data[0],psdCumm = psd_1_data[1],
                  distributeMass = False)

# Define two type of engines 
if contact_type == 0:
    print 'No particle rolling!'
    print '\n'
    sp0.toSimulation(material = particles_mat_0)
    O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_FrictPhys()],
            [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        consolidation,
        NewtonIntegrator(damping=.2)
    ]
    
if contact_type == 1:
    print 'Relative particle rolling is considered!'
    print '\n'
    sp0.toSimulation(material = particles_mat_1)
    O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
        InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
            [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
                useIncrementalForm=True, 
                always_use_moment_law=True,  
                label='cohesiveLaw')]
        ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        consolidation,
        NewtonIntegrator(damping=.2),
    ]

# Running the task
O.dt = utils.PWaveTimeStep()
O.usesTimeStepper = True

for i in xrange(10):
    O.run(200, True)
    radius_list = []
    for b in O.bodies:
        if isinstance(b.shape,Sphere):
            init_particles_rad = b.shape.radius
            radius_list.append(init_particles_rad)
    min_radius = min(radius_list)
    max_radius = max(radius_list)
    print 'The initial min radius is: ',min_radius
    print 'The initial max radius is: ',max_radius
    print '\n'
O.wait()

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.