New question #706797 on Yade: https://answers.launchpad.net/yade/+question/706797
Hello, I create a sphere packing position in the 1st script. However in second script, I import the position from the 1st script. I reduce the sphere value. I am confuse on why once the sphere hit the funnel the simulation stop right away. ### 1st script ### import random import math from yade import geom, pack, utils, plot, ymport, export import numpy as np # Define cylinder with funnel parameters center = (0, 0, 0) diameter = 0.102 height = 0.18 # create cylindrical body with radius 0.102 m and height 0.064 m cylinder = geom.facetCylinder(center=center, radius=diameter/2, height=height, segmentsNumber=80, wallMask=6) # add cylinder to simulation O.bodies.append(cylinder) # Define cylinder with funnel parameters center1 = (0,0,height/2) dBunker = 0.4 dOutput = 0.102 hBunker = 0 hOutput = 0.15 hPipe = 0 # create funnel as a bunker with diameter 0.102 m, height 0.064 m funnel = geom.facetBunker(center=center1, dBunker=dBunker, dOutput=dOutput, hBunker=hBunker,hOutput=hOutput, hPipe=hPipe, segmentsNumber=80, wallMask=4) # add funnel to simulation O.bodies.append(funnel) def makeRandomBlocks(diameter, gridsize, minz, numblocks): bunkerLimit = diameter/2 # # Make grid blocks inside a bunker # grid = np.arange(-bunkerLimit, bunkerLimit, gridsize) inGrid = [] for yi in grid: for xi in grid: c = [[xi, yi], [xi+gridsize, yi], [xi+gridsize, yi+gridsize], [xi, yi+gridsize]] # # Check if the block corners are outside of the bunker or not # out = False for p in c: r = (p[0]**2 + p[1]**2)**0.5 if r > bunkerLimit: out = True # # If the block is inside the bunker, keep it # if not out: inGrid.append(c) layer = math.ceil(numblocks / len(inGrid)) blocks = [] for l in range(layer): for g in inGrid: zi = minz + l*gridsize p1 = g[0].copy() p1.append(zi) p2 = g[2].copy() p2.append(zi+gridsize) blocks.append([p1, p2]) random.shuffle(blocks) return blocks minZ = 0.35 dbunker = 0.4 gridSize = dbunker/8 numblocks = 51 blockList = makeRandomBlocks(dbunker, gridSize, minZ, numblocks) for i in range(numblocks): print("Making cloud block", i+1, "/", numblocks) corner = blockList.pop() sp = pack.SpherePack() # 15.75 mm 13 particles if (i < 13): n1 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.01575/2, num=1) # 11 mm 51 particles n2 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.011/2, num=1) # 7.125 mm 51x11 = 561 particles n3 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.007125/2, num=11) # 3.555 mm 51x100 = 5,100 particles n4 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.003555/2, num=100) # 1.60 mm 51x1400 = 71,400 particles n5 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.0016/2, num=1400) # 0.8 mm 51x4140 = 211,140 particles n6 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.0008/2, num=4140) sp.toSimulation() export.text("testCloud.txt") ################################################## ### 2nd script ### import random import math from yade import geom, pack, utils, plot, ymport, export import numpy as np # Define cylinder with funnel parameters center = (0, 0, 0) diameter = 0.102 height = 0.18 # create cylindrical body with radius 0.102 m and height 0.064 m cylinder = geom.facetCylinder(center=center, radius=diameter/2, height=height, segmentsNumber=80, wallMask=6) # add cylinder to simulation O.bodies.append(cylinder) # Define cylinder with funnel parameters center1 = (0,0,height/2) dBunker = 0.4 dOutput = 0.102 hBunker = 0 hOutput = 0.15 hPipe = 0 # create funnel as a bunker with diameter 0.102 m, height 0.064 m funnel = geom.facetBunker(center=center1, dBunker=dBunker, dOutput=dOutput, hBunker=hBunker,hOutput=hOutput, hPipe=hPipe, segmentsNumber=80, wallMask=4) # add funnel to simulation O.bodies.append(funnel) # add sphere packing O.bodies.append(ymport.textExt('testCloud.txt',format='x_y_z_r')) # materials Properties gravel = CohFrictMat(young = 1e7, poisson = 0.25, density = 27000, label = 'gravel') asphalt_binder = CohFrictMat(young = 1e7, poisson = 0.25, density = 10600, frictionAngle = radians(40), normalCohesion = 2e5, shearCohesion = 2e5, label = 'asphalt_binder') # add properties O.materials.append(gravel) O.materials.append(asphalt_binder) # give color and properties to shpere for body in O.bodies: if not isinstance(body.shape, Sphere): continue if body.shape.radius == 0.01575/2 : body.shape.color = (0,0,1) #blue body.material = gravel if body.shape.radius == 0.011/2: body.shape.color = (1,0,0) #red body.material = gravel if body.shape.radius == 0.007125/2: body.shape.color = (0,1,0) #green body.material = gravel if body.shape.radius == 0.003555/2: body.shape.color = (1,1,0) #yellow body.material = gravel if body.shape.radius == 0.00160/2 : body.shape.color = (1,0,1) #magenta body.material = gravel if body.shape.radius == 0.0008/2 : body.shape.color = (0,0,0) #black body.material = asphalt_binder O.engines = [ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]), InteractionLoop( # handle sphere+sphere and facet+sphere collisions [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()] ), NewtonIntegrator(gravity=(0, 0, -981), damping=0.4), # call the checkUnbalanced function (defined below) every 2 seconds PyRunner(command='checkUnbalanced()', realPeriod=2), # call the addPlotData function every 200 steps PyRunner(command='addPlotData()', iterPeriod=100) ] O.dt = 0.1 * PWaveTimeStep() # enable energy tracking; any simulation parts supporting it # can create and update arbitrary energy types, which can be # accessed as O.energy['energyName'] subsequently O.trackEnergy = True # if the unbalanced forces goes below .05, the packing # is considered stabilized, therefore we stop collected # data history and stop def checkUnbalanced(): if unbalancedForce() < 0.01: O.pause() export.text("gra_depo_pos.txt") plot.saveDataTxt('bbb.txt') # plot.saveGnuplot('bbb') is also possible # 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(), **O.energy) # define how to plot data: 'i' (step number) on the x-axis, unbalanced force # on the left y-axis, all energies on the right y-axis # (O.energy.keys is function which will be called to get all defined energies) # None separates left and right y-axis plot.plots = {'i': ('unbalanced', None, O.energy.keys)} # show the plot on the screen, and update while the simulation runs plot.plot() O.saveTmp() ################################################## -- 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