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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp