yade-users team mailing list archive

[Ruidong LI <question707183@xxxxxxxxxxxxxxxxxxxxx>]

New question #707183 on Yade:
https://answers.launchpad.net/yade/+question/707183

Hi! I want to model the cylindrical compression and follow the procedures as follows:
1.	Import geometry from microCT
2.	Create clumps using the Euclidean 3D method mentioned in the 'Clump' code [1]
3.	Equilibrate within rigid boundaries until unbalanced force is zero
4.	Apply gravity and leave run until unbalanced force is zero
5.	Check void ratio.

But the problems are so werid and summarized as follows:
1. Sample height decrease
2. High unbalanced force ratio after applying gravity.

My MWE is as follows:

from __future__ import print_function
from yade import pack, ymport, qt
import gts, os.path, locale, plot, random
import numpy as np
import math

locale.setlocale(
        locale.LC_ALL, 'en_US.UTF-8'
)  #gts is locale-dependend.  If, for example, german locale is used, gts.read()-function does not import floats normally

###

############################################
###                  1                   ###
###   DEFINING VARIABLES AND MATERIALS   ###
############################################

# The following lines are used parameter definitions
readParamsFromTable(
	# directory	
	dir_clump = '/home/kyle/Downloads/yadetxtClumps20_1.16.txt', # positions and radius of spherical particles in the clump
	dir_vtk = '/home/kyle/Yade/install/bin/lentils/vtk/20', # directory of storing vtk files
	dir_txt = '/home/kyle/Yade/install/bin/lentils/txt/20/Num20_initialpacking.txt', # directory of storing txt files
	dir_xml = '/home/kyle/Yade/install/bin/lentils/xml/20', # directory of storing xml files
        
        # material parameters
        young_w = 5e9, # Young's modulus of plates and walls
        young_g = 1.8e6, # Young's modulus of grains [1]
        den_w = 7874, # density of plates and walls
        den_g = 980, # density of grains 
        poi_w = 0.3, # possion's ratio of plates and walls
        poi_g = 0.25, # possion's ratio of grains
        friAngle_w = 0.5,#radians(15), # friction angle of plates and walls 
        friAngle_g = 0.5,#radians(29), # friction angle of grains

        # Parameters of cylindrical walls
        x_cyl = 0.0547, # x-coordinate of center of cylindrical walls
        y_cyl = 0.0535, # y-coordinate of center of cylindrical walls
        z_cyl = 0, # z-coordinate of center of cylindrical walls
        r_cyl = 0.0358, # radius of the cylinder
        h_cyl = 0.14, # height of the cylinder
        n_cyl = 100, # number of boxes forming cylinder
        t_cyl = 0.0001, # thickness of the boxes
        
        #Parameters of lower and upper plates (square box)
        w_plate = 0.05, # width of plates
        t_plate = 0.0001, # thickness of plates
        
)
from yade.params.table import *

# create materials for spheres and walls
wallMat = O.materials.append(FrictMat(young = young_w, poisson = poi_w, frictionAngle = friAngle_w, density = den_w, label = 'walls'))
sphereMat = O.materials.append(FrictMat(young = young_g, poisson = poi_g, frictionAngle = friAngle_g, density = den_g, label = 'spheres'))

###############################################
###                    2                    ###
###   IMPORTING GRAINS AND CREATING CLUMPS  ###
###############################################

# import clumps
spheres = ymport.textClumps(dir_clump,material=sphereMat,discretization=0)

# assign unique color for each clump
currentClumpId = O.bodies[0].clumpId
color = randomColor()
for b in O.bodies:
    if b.clumpId != currentClumpId:#change color each time clumpId changes
        currentClumpId = b.clumpId
        color = randomColor()
    if isinstance(b.shape,Sphere):#colorize spheres
        b.shape.color = color

# count the number of clumps 
bodyList = []
for i in O.bodies:
 if i.isClump:
  bodyList.append(i.id)
print(len(bodyList))

###################################################################
###                              3                              ###
###   CREATING CYLINDRICAL WALLS, BOTTOM PLATE, AND TOP PLATE   ###
###################################################################

h = max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])
#top_plate = O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,h),radius=r_cyl*1.2,height=0,segmentsNumber=40,color=(1,1,1)))
#bottom_plate = O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,0),radius=r_cyl*1.2,height=0,segmentsNumber=40,color=(1,1,1)))
cyl_wallRigid = O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,h/2),radius=r_cyl,height=h,segmentsNumber=100,color=(0,0,1)))

# calculate void ratio
def myClumpVoidRatio(VolTot, density):
	mass=sum([ b.state.mass for b in O.bodies if b.isClump])
	VolSolid = mass/density
	poro = (VolTot-VolSolid)/VolSolid
	return poro	
h_cyl = max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])
Vol = math.pi * h_cyl * r_cyl * r_cyl
poro_ini = myClumpVoidRatio(Vol, den_g)
print('Initial void ratio:', str(poro_ini))

############################
###           4          ###
###   DEFINING ENGINES   ###
############################

O.engines = [
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
	InteractionLoop([Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
	PyRunner(command='cycleClumps()', iterPeriod = 100, label='cycle'),
	PyRunner(command='addPlotData()', iterPeriod=100),
	VTKRecorder(iterPeriod = 1000, recorders = ['all'], fileName = dir_vtk + '/'),
	NewtonIntegrator(gravity = [0, 0, 0], damping = 0.4, label="NI")
]
O.dt= 0.5 * PWaveTimeStep()

# cycle clumps until the unbalance force ratio is sufficient small
global n_ini
n_ini = 0

def myClumpPorosity(VolTot, density):
  mass=sum([ b.state.mass for b in O.bodies if b.isClump])
  return (VolTot -mass/density)/VolTot

def cycleClumps():
	h_cyl = max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])
	Vol = math.pi * h_cyl * r_cyl * r_cyl
	poro_ini = myClumpVoidRatio(Vol, den_g)
	global n_ini
	print(O.iter, unbalancedForce(), n_ini, max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]), poro_ini)
	f = open(dir_txt, 'a')
	f.write('{} {} {} {} \n'.format(O.iter, unbalancedForce(), n_ini, max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])))
	f.close()
	if unbalancedForce() < 1e-3:
		if n_ini == 0:
			NI.gravity = Vector3(0,0,-9.81)
			n_ini = n_ini + 1
	if n_ini == 1:
		if O.iter > 5000000:
		#if unbalancedForce() > 1e-2:
			O.pause()
			O.save('firstcalculation.xml')

# collect history of data which will be plotted
def addPlotData():
	# each item is given a names, by which it can be the unsed in plot.plots
	# the **O.energy converts dictionary-like O.energy to plot.addData arguments
	plot.addData(i=O.iter, unbalanced=unbalancedForce(), coordNum=avgNumInteractions(considerClumps=True), **O.energy)

##########################
###          5         ###
###   RUN SIMULATION   ###
##########################
qt.View()
rr = yade.qt.Renderer()
rr.shape = True
rr.intrPhys = True

plot.plots = {'i': ('unbalanced', None, 'coordNum')}
plot.plot()	
O.run()

[1]https://github.com/ElsevierSoftwareX/SOFTX-D-21-00045

The initial txt named 'yadetxtClumps20_1.16.txt' can be downloaded from: https://drive.google.com/file/d/17yyABYuXUdMPnr2x-QqQwhA_ParjnHMw/view?usp=sharing


-- 
You received this question notification because your team yade-users is
an answer contact for Yade.