yade-users team mailing list archive

[Fu zuoguang <question258445@xxxxxxxxxxxxxxxxxxxxx>]

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

Fu zuoguang gave more information on the question:
codes:
#########################################################################################
# unicode: UTF-8
Filename='3D-triaxial compression test with a given PSD'
from yade import pack,os,utils
import math

'''
The purpose of this simulation is to run sucessfully a 2D-biaxial compression
with a pre-define PSD in Yade DEM platform and then to compare its results to
lab test achievements.

The whole processing of this simulation can be divided into 3 parts, which can
be described as follow:

[1].material defination. Of course, there are just two types of materials needed
to be defined, the first one is the particles-material, cohefrict-cohefrict material
(incohesive but can consider rolling mechanism)
and the other is the walls-material, frict-frict material(incohesive).

[2].Bodies generation. In this simulation, just two types of bodies should be considered,
particles which here are simplified as idealistic spheres and walls which are introduced
for representing the sample's boundary. These two here are treated as rigid body. And the
particles are generated with a given PSD.

[3].Running the simulation. The whole process of this simulation composes of just two parts,
which can be called, of course, the isotropic consolidation process and the triaxial
compression process. So this step is to define the calculation engines for consolidation.
Despite that the algorithms of contact detection in Yade (it can be said that in all of
the DEM software here) and of random packing in certain space make it difficult to satisfy
that the initial porosity of a given simulation sample at the start of the consolidtaion
is as the same as it of one in lab, some mathematics methods can be employed as a approximate
solution, such as radius expansion menthod, to obtain a appropriate porosity with a certain
confine pressure.Some more details are in the below section.
'''

##################################  materials defination ###############################
'''
O.materials.append(FrictMat(young=1e9,poisson=0.5,frictionAngle=radians(20),
                    density=2600,label='particles_mat'))
O.materials.append(FrictMat(young=1e9,poisson=0.5,frictionAngle=0,
                    density=0,label='walls_mat'))
'''
# Read the simulation condition.
condition_data = []
f = file('condition.txt', 'r')
for line in f.readlines():
    condition_data.append(line.strip('\n').split(':'))
for i in xrange(len(condition_data)):    
    if i < 2 :
        condition_data[i][1] = condition_data[i][1].split(',')
        for j in xrange(len(condition_data[i][1])):
            condition_data[i][1][j] = float(condition_data[i][1][j])
        condition_data[i][1] = tuple(condition_data[i][1])
    elif i == 2 :
        condition_data[i][1] = int(condition_data[i][1])
    elif i == 3 :
        condition_data[i][1] = bool(float(condition_data[i][1]))
    else:
        condition_data[i][1] = float(condition_data[i][1])
f.close()
# Particles-material defination
material_list_1 = []
material_list_1.append(8e8)             # list[0]='young module'
material_list_1.append(0.5)             # list[1]='poisson ratio'
material_list_1.append(0.59)            # list[2]='friction Angle'
material_list_1.append(2650)            # list[3]='density'
material_list_1.append(0.3)             # list[4]='rolling coefficient'
material_list_1.append(1.0)             # list[5]='considering plastic moment of rolling'
compFricDegree = material_list_1[2]     # for radius expansion
particles_mat = O.materials.append(CohFrictMat( young = material_list_1[0],
                                                poisson = material_list_1[1],
                                                frictionAngle = material_list_1[2],
                                                density  =material_list_1[3],
                                                isCohesive = False,
                                                momentRotationLaw = True,
                                                alphaKr = material_list_1[4],
                                                etaRoll = material_list_1[5]))

# walls-material defination
material_list_2 = []
material_list_2.append(1e9)            # list[0]='young module'
material_list_2.append(0.5)             # list[1]='poisson ratio'
material_list_2.append(radians(0))      # list[2]='friction Angle'
material_list_2.append(0)               # list[3]='density'
walls_mat = O.materials.append(FrictMat(young=material_list_2[0],
                                        poisson=material_list_2[1],
                                        frictionAngle=material_list_2[2],
                                        density=material_list_2[3]))

##################################  bodies defination
###############################

# walls generation.
mincorner = condition_data[0][1]
maxcorner = condition_data[1][1]
walls_width = 1e-9
# In 2D simulation, spheres generation corners are:
# genmincorner = [mincorner[0],mincorner[1],(maxcorner[2]-mincorner[2])/2.0]
# genmaxcorner = [maxcorner[0],maxcorner[1],(maxcorner[2]-mincorner[2])/2.0]
# In 3D simulation, spheres generation corners are:
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)
left_wall_posi = O.bodies[0].state.pos[0]        # pos[0]--x direction
right_wall_posi = O.bodies[1].state.pos[0]
bottom_wall_posi = O.bodies[2].state.pos[1]      # pos[1]--y direction
top_wall_posi = O.bodies[3].state.pos[1]
back_wall_posi = O.bodies[4].state.pos[2]        # pos[2]--z direction
front_wall_posi = O.bodies[5].state.pos[2]
x_length = abs(right_wall_posi - left_wall_posi)
y_length = abs(top_wall_posi - bottom_wall_posi)
z_length = abs(front_wall_posi - back_wall_posi)
# In 3D simulation:
sample_volume = x_length * y_length * z_length
# In 2D simualtion:
# sample_volume = x_length * y_length
print 'The x_length of sample is: ',x_length
print 'The y_length of sample is: ',y_length
print 'The z_length of sample is: ',z_length
print 'The volume of sample is: ',sample_volume,'\n'

# particles generation with a given PSD curve. The objective is to ensure
# that all the particles generated in packbox can have different unique
# sizes and all these sizes can fix approximately the PSD data of realistic
# granular materials. The programs for this purpose are composed of 7 parts,
# which are as follows:

# The first step: read PSD data for pre-generation.
f = file('psd.txt','r')
read_PSD_data_str = []; read_PSD_data_flo = []
for line in f.readlines():
    datalist = line.split()
    read_PSD_data_str.append(datalist)
for i in xrange(1,len(read_PSD_data_str),1):
    read_PSD_data_flo.append([])
    for j in read_PSD_data_str[i]:
        j = float(j)
        read_PSD_data_flo[i-1].append(j)

# The second step: particles generating in makecloud with this given PSD.
PSD_diameters_list = []; PSD_frequencies_list = [];PSD = []
for i in xrange(len(read_PSD_data_flo)):
    PSD_diameters_list.append(read_PSD_data_flo[i][0])
    PSD_frequencies_list.append(read_PSD_data_flo[i][1])
PSD.append(genmincorner)
PSD.append(genmaxcorner)
PSD.append(condition_data[2][1])
PSD.append(PSD_diameters_list)
PSD.append(PSD_frequencies_list)
sp0 = pack.SpherePack();
sp0.makeCloud(PSD[0],PSD[1],num = PSD[2],psdSizes = PSD[3],
              psdCumm = PSD[4],distributeMass=True)

# The third step: particles generating from makecloud to yade.
sp0.toSimulation(material=particles_mat)
radius_list = []
for b in O.bodies:
    if isinstance(b.shape,Sphere):
        init_particles_rad = b.shape.radius
        radius_list.append(init_particles_rad)
'''
# The fourth step: blockedDOFs for 2D simulation tasks.
for b in O.bodies:
    if isinstance(b.shape,Sphere):
        b.state.blockedDOFs='zXY'
'''
# Making an output to test the particles number and the min or max radii.
min_diameter = 2*min(radius_list);
max_diameter = 2*max(radius_list)
print 'The total number of particles in simulation is: ',len(radius_list)
print 'The initial max diameter is: ',max_diameter
print 'The initial min diameter is: ',min_diameter,'\n'

# Making post-processing for PSD, which may include:
# [1]. particle diameter VS relative frequency curve by munber accumulation.
# [2]. particle diameter VS cumulative frequency curve by munber accumulation.
# [3]. particle diameter VS cumulative frequency curve by volume(mass) accumulation.
particle_dia_VS_num_rel = []
particle_dia_VS_num_cum = []
particle_dia_VS_mass_cum = []
maxd = max_diameter
mind = min_diameter
mund = len(PSD_diameters_list)-1
for i in xrange(mund):
    particle_dia_VS_num_rel.append([])
    particle_dia_VS_num_cum.append([])
    particle_dia_VS_mass_cum.append([])
    lower_level_limit = mind + (maxd-mind)*i/mund
    upper_level_limit = mind + (maxd-mind)*(i+1)/mund
    for j in radius_list:
        if lower_level_limit <= 2*j and 2*j <= upper_level_limit:
            particle_dia_VS_num_rel[i].append(2*j)
        if 2*j <= upper_level_limit:
            particle_dia_VS_num_cum[i].append(2*j)
            particle_dia_VS_mass_cum[i].append((4/3.0)*(math.pi)*j*j*j)

# Making output of PSD data, which may include:
# [1]. the data of diameter VS relative frequency curve in number based.
# [2]. the data of diameter VS commulative frequency curve in number based.
# [3]. mean radii of particles according to cummulative number.
# [4]. the data of diameter VS commulative frequency curve in mass based.
# [5]. mean radii of particles according to cummulative mass.
# [6]. the porosity of this sample.
print 'the data of diameter VS relative frequency curve in number based are: '
for i in xrange(mund):
    print float(len(particle_dia_VS_num_rel[i]))/len(radius_list)
print '\n'

print 'the data of diameter VS commulative frequency curve in number based are: '
for i in xrange(mund):
    print float(len(particle_dia_VS_num_cum[i]))/len(radius_list)
print '\n'

print 'mean radii of particles according to cummulative number is: '
mean_radii_size = sum(radius_list)/len(radius_list)
print mean_radii_size,'\n'

ori_particles_sum_volume = []
for i in xrange(len(radius_list)):
    orirad = radius_list[i]
    ori_particle_volume = ((4/3.0)*(math.pi)*(orirad**3))
    ori_particles_sum_volume.append(ori_particle_volume)
initial_particles_vol = float(sum(ori_particles_sum_volume))
print 'the data of diameter VS commulative frequency curve in mass based are: '
for i in xrange(mund):
    print float(sum(particle_dia_VS_mass_cum[i]))/initial_particles_vol
print '\n'

print 'mean radii of particles according to cummulative mass is: '
mean_radius_vol = ((initial_particles_vol)/((4/3.0)*(math.pi)*len(radius_list)))**(1.0/3.0)
print mean_radius_vol,'\n'

print 'the porosity of this sample is: '
initial_porosity = 1-(initial_particles_vol)/sample_volume
print initial_porosity,'\n'

# Making an output for getting all the desired data,which include:
# [1]. PSD-data.dat
# [2]. radius-data.dat
f=file("PSD-information.csv",'a')
f.write('%-12s,%-12s,%-12s,%-12s,\n'%('Diameter','Frequency',
                                    'Cum-number','Cum-mass'))
for i in xrange(mund):
    lower_level_limit = mind + (maxd-mind)*i/mund
    upper_level_limit = mind + (maxd-mind)*(i+1)/mund
    a = ( lower_level_limit + upper_level_limit )/2
    b = float(len(particle_dia_VS_num_rel[i]))/len(radius_list)
    c = float(len(particle_dia_VS_num_cum[i]))/len(radius_list)
    d = float(sum(particle_dia_VS_mass_cum[i]))/initial_particles_vol
    f.write('%12.5E,%12.5E,%12.5E,%12.5E,\n'%(a,b,c,d))
f.close()

############################## consolidation engines defination #####################
sample_density = condition_data[3][1]
confine_pressure = condition_data[4][1]
strain_rate = condition_data[5][1]
final_axial_strain = condition_data[6][1]
print '#####################################################'

consolidation = TriaxialStressController(
# defining the ids of walls, of which the normal directions are respectively parallel to x,y,z-coordinate.
    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,
# In consolidation process, choosing the 'radius expansion' method which can provide the chance that
# all the radius of particles may increase with the certain Multiplier until the target confine-pressure
# or porosity are attained. At this time, we choose the confine-pressure as the target.
    maxMultiplier = 1.01,           # spheres growing factor (fast growth)
    finalMaxMultiplier = 1.001,     # spheres growing factor (slow growth)
    thickness=walls_width,
    internalCompaction=True,        # 'False' here may close the task

# switch stress/strain control using a bitmask. What is a bitmask ?
# Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
# Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
# to put it differently, the mask is the integer whose binary representation is xyz, i.e.
# "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
    stressMask=7,
        goal1 = - confine_pressure,    # confine-pressure applied in x-derection.
        goal2 = - confine_pressure,    # confine-pressure applied in y-derection.
        goal3 = - confine_pressure,    # confine-pressure applied in z-derection.
# Defining the max velocity of walls to make sure the isotropic conpression in the whole process
#       max_vel = 10,
)

# Define consolidation results output
f=file("consolidation.csv",'a')
f.write('%-12s,%-12s,%-12s,\n'%('Iter','confine-prs','porosity'))
f.close()

def consolidation_results():
# Define mean sigma as the confinment pressure.
    sigma_1 = -consolidation.stress(consolidation.wall_top_id)[1]
    sigma_2 = -consolidation.stress(consolidation.wall_right_id)[0]
    sigma_3 = -consolidation.stress(consolidation.wall_front_id)[2]
    meanS = (sigma_1 + sigma_2 + sigma_3)/3.0
# Define variable porosity in the process of consolidation.
    sum_particles_volume = []; max_particles_rad = []
    for b in O.bodies:
        if isinstance(b.shape,Sphere):
            rad = b.shape.radius
            single_particles_volume = (4.0/3.0)*(math.pi)*(rad**3)
            sum_particles_volume.append(single_particles_volume)
            max_particles_rad.append(rad)
    total_particles_volume = sum(sum_particles_volume)
    max_particles_radii = max(max_particles_rad)
    print 'the max radii of particles is: ',max_particles_radii
    print 'the total particle volume is: ',total_particles_volume
    porosity = 1 - (total_particles_volume/sample_volume)
# Make all the things above outputing
    f=file("consolidation.csv",'a')
    f.write('%12.5E,%12.5E,%12.5E,\n'%(O.iter,meanS,porosity)),
    f.close()

# Define simulation engine
O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
    InteractionLoop(
# Box-sphere interactions will be the simple normal-shear law, we use ScGeom for them
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
# Boxes will be frictional (FrictMat), so the sphere-box physics is FrictMat vs. CohFrictMat.
        [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
# Finally, two different contact laws for sphere-box and sphere-sphere
        [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),
#   PyRunner(command='consolidation_results()',iterPeriod=200)
]

# The class 'PWaveTimeStep' is used as 'dt'.
O.dt = utils.PWaveTimeStep()
O.usesTimeStepper = True

for i in xrange(10):
    O.run(100,True)
    consolidation_results()
O.wait()

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