Examples with tutorial¶
The online version of this tutorial contains embedded videos.
Bouncing sphere¶
Following example is in file doc/sphinx/tutorial/01-bouncing-sphere.py.
# basic simulation showing sphere falling ball gravity,
# bouncing against another sphere representing the support
# DATA COMPONENTS
# add 2 particles to the simulation
# they the default material (utils.defaultMat)
O.bodies.append(
[
# fixed: particle's position in space will not change (support)
sphere(center=(0, 0, 0), radius=.5, fixed=True),
# this particles is free, subject to dynamics
sphere((0, 0, 2), .5)
]
)
# FUNCTIONAL COMPONENTS
# simulation loop -- see presentation for the explanation
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()], # collision geometry
[Ip2_FrictMat_FrictMat_FrictPhys()], # collision "physics"
[Law2_ScGeom_FrictPhys_CundallStrack()] # contact law -- apply forces
),
# Apply gravity force to particles. damping: numerical dissipation of energy.
NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.1)
]
# set timestep to a fraction of the critical timestep
# the fraction is very small, so that the simulation is not too fast
# and the motion can be observed
O.dt = .5e-4 * PWaveTimeStep()
# save the simulation, so that it can be reloaded later, for experimentation
O.saveTmp()
Gravity deposition¶
Following example is in file doc/sphinx/tutorial/02-gravity-deposition.py.
# gravity deposition in box, showing how to plot and save history of data,
# and how to control the simulation while it is running by calling
# python functions from within the simulation loop
# import yade modules that we will use below
from yade import pack, plot
# create rectangular box from facets
O.bodies.append(geom.facetBox((.5, .5, .5), (.5, .5, .5), wallMask=31))
# create empty sphere packing
# sphere packing is not equivalent to particles in simulation, it contains only the pure geometry
sp = pack.SpherePack()
# generate randomly spheres with uniform radius distribution
sp.makeCloud((0, 0, 0), (1, 1, 1), rMean=.05, rRelFuzz=.5)
# add the sphere pack to the simulation
sp.toSimulation()
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, -9.81), 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 = .5 * 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() < .05:
O.pause()
plot.saveDataTxt('bbb.txt.bz2')
# 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()
Oedometric test¶
Following example is in file doc/sphinx/tutorial/03-oedometric-test.py.
# gravity deposition, continuing with oedometric test after stabilization
# shows also how to run parametric studies with yade-batch
# The components of the batch are:
# 1. table with parameters, one set of parameters per line (ccc.table)
# 2. readParamsFromTable which reads respective line from the parameter file
# 3. the simulation muse be run using yade-batch, not yade
#
# $ yade-batch --job-threads=1 03-oedometric-test.table 03-oedometric-test.py
#
# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(rMean=.05, rRelFuzz=.3, maxLoad=1e6, minLoad=1e4)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *
# create box with free top, and ceate loose packing inside the box
from yade import pack, plot
O.bodies.append(geom.facetBox((.5, .5, .5), (.5, .5, .5), wallMask=31))
sp = pack.SpherePack()
sp.makeCloud((0, 0, 0), (1, 1, 1), rMean=rMean, rRelFuzz=rRelFuzz)
sp.toSimulation()
O.engines = [
ForceResetter(),
# sphere, facet, wall
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(), Bo1_Wall_Aabb()]),
InteractionLoop(
# the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.5),
# the label creates an automatic variable referring to this engine
# we use it below to change its attributes from the functions called
PyRunner(command='checkUnbalanced()', realPeriod=2, label='checker'),
]
O.dt = .5 * PWaveTimeStep()
# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time
# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
# at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
if O.iter < 5000:
return
# the rest will be run only if unbalanced is < .1 (stabilized packing)
if unbalancedForce() > .1:
return
# add plate at the position on the top of the packing
# the maximum finds the z-coordinate of the top of the topmost particle
O.bodies.append(wall(max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]), axis=2, sense=-1))
global plate # without this line, the plate variable would only exist inside this function
plate = O.bodies[-1] # the last particles is the plate
# Wall objects are "fixed" by default, i.e. not subject to forces
# prescribing a velocity will therefore make it move at constant velocity (downwards)
plate.state.vel = (0, 0, -.1)
# start plotting the data now, it was not interesting before
O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=200)]
# next time, do not call this function anymore, but the next one (unloadPlate) instead
checker.command = 'unloadPlate()'
def unloadPlate():
# if the force on plate exceeds maximum load, start unloading
if abs(O.forces.f(plate.id)[2]) > maxLoad:
plate.state.vel *= -1
# next time, do not call this function anymore, but the next one (stopUnloading) instead
checker.command = 'stopUnloading()'
def stopUnloading():
if abs(O.forces.f(plate.id)[2]) < minLoad:
# O.tags can be used to retrieve unique identifiers of the simulation
# if running in batch, subsequent simulation would overwrite each other's output files otherwise
# d (or description) is simulation description (composed of parameter values)
# while the id is composed of time and process number
plot.saveDataTxt(O.tags['d.id'] + '.txt')
O.pause()
def addPlotData():
if not isinstance(O.bodies[-1].shape, Wall):
plot.addData()
return
Fz = O.forces.f(plate.id)[2]
plot.addData(Fz=Fz, w=plate.state.pos[2] - plate.state.refPos[2], unbalanced=unbalancedForce(), i=O.iter)
# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots = {'i': ('unbalanced',), 'w': ('Fz',)}
plot.plot()
O.run()
# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()
Batch table¶
To run the same script doc/sphinx/tutorial/03-oedometric-test.py in batch mode to test different parameters, execute command yade-batch 03-oedometric-test.table 03-oedometric-test.py, also visit page http://localhost:9080 to see the batch simulation progress.
rMean rRelFuzz maxLoad
.05 .1 1e6
.05 .2 1e6
.05 .3 1e6
Periodic simple shear¶
Following example is in file doc/sphinx/tutorial/04-periodic-simple-shear.py.
# encoding: utf-8
# script for periodic simple shear test, with periodic boundary
# first compresses to attain some isotropic stress (checkStress),
# then loads in shear (checkDistorsion)
#
# the initial packing is either regular (hexagonal), with empty bands along the boundary,
# or periodic random cloud of spheres
#
# material friction angle is initially set to zero, so that the resulting packing is dense
# (sphere rearrangement is easier if there is no friction)
#
# setup the periodic boundary
O.periodic = True
O.cell.hSize = Matrix3(2, 0, 0, 0, 2, 0, 0, 0, 2)
from yade import pack, plot
# the "if 0:" block will be never executed, therefore the "else:" block will be
# to use cloud instead of regular packing, change to "if 1:" or something similar
if 0:
# create cloud of spheres and insert them into the simulation
# we give corners, mean radius, radius variation
sp = pack.SpherePack()
sp.makeCloud((0, 0, 0), (2, 2, 2), rMean=.1, rRelFuzz=.6, periodic=True)
# insert the packing into the simulation
sp.toSimulation(color=(0, 0, 1)) # pure blue
else:
# in this case, add dense packing
O.bodies.append(pack.regularHexa(pack.inAlignedBox((0, 0, 0), (2, 2, 2)), radius=.1, gap=0, color=(0, 0, 1)))
# create "dense" packing by setting friction to zero initially
O.materials[0].frictionAngle = 0
# simulation loop (will be run at every step)
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
# interaction loop
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(damping=.4),
# run checkStress function (defined below) every second
# the label is arbitrary, and is used later to refer to this engine
PyRunner(command='checkStress()', realPeriod=1, label='checker'),
# record data for plotting every 100 steps; addData function is defined below
PyRunner(command='addData()', iterPeriod=100)
]
# set the integration timestep to be 1/2 of the "critical" timestep
O.dt = .5 * PWaveTimeStep()
# prescribe isotropic normal deformation (constant strain rate)
# of the periodic cell
O.cell.velGrad = Matrix3(-.1, 0, 0, 0, -.1, 0, 0, 0, -.1)
# when to stop the isotropic compression (used inside checkStress)
limitMeanStress = -5e5
# called every second by the PyRunner engine
def checkStress():
# stress tensor as the sum of normal and shear contributions
# Matrix3.Zero is the intial value for sum(...)
stress = getStress().trace() / 3.
print('mean stress', stress)
# if mean stress is below (bigger in absolute value) limitMeanStress, start shearing
if stress < limitMeanStress:
# apply constant-rate distorsion on the periodic cell
O.cell.velGrad = Matrix3(0, 0, .1, 0, 0, 0, 0, 0, 0)
# change the function called by the checker engine
# (checkStress will not be called anymore)
checker.command = 'checkDistorsion()'
# block rotations of particles to increase tanPhi, if desired
# disabled by default
if 0:
for b in O.bodies:
# block X,Y,Z rotations, translations are free
b.state.blockedDOFs = 'XYZ'
# stop rotations if any, as blockedDOFs block accelerations really
b.state.angVel = (0, 0, 0)
# set friction angle back to non-zero value
# tangensOfFrictionAngle is computed by the Ip2_* functor from material
# for future contacts change material (there is only one material for all particles)
O.materials[0].frictionAngle = .5 # radians
# for existing contacts, set contact friction directly
for i in O.interactions:
i.phys.tangensOfFrictionAngle = tan(.5)
# called from the 'checker' engine periodically, during the shear phase
def checkDistorsion():
# if the distorsion value is >.3, exit; otherwise do nothing
if abs(O.cell.trsf[0, 2]) > .5:
# save data from addData(...) before exiting into file
# use O.tags['id'] to distinguish individual runs of the same simulation
plot.saveDataTxt(O.tags['id'] + '.txt')
# exit the program
#import sys
#sys.exit(0) # no error (0)
O.pause()
# called periodically to store data history
def addData():
# get the stress tensor (as 3x3 matrix)
stress = sum(normalShearStressTensors(), Matrix3.Zero)
# give names to values we are interested in and save them
plot.addData(exz=O.cell.trsf[0, 2], szz=stress[2, 2], sxz=stress[0, 2], tanPhi=(stress[0, 2] / stress[2, 2]) if stress[2, 2] != 0 else 0, i=O.iter)
# color particles based on rotation amount
for b in O.bodies:
# rot() gives rotation vector between reference and current position
b.shape.color = scalarOnColorScale(b.state.rot().norm(), 0, pi / 2.)
# define what to plot (3 plots in total)
## exz(i), [left y axis, separate by None:] szz(i), sxz(i)
## szz(exz), sxz(exz)
## tanPhi(i)
# note the space in 'i ' so that it does not overwrite the 'i' entry
plot.plots = {'i': ('exz', None, 'szz', 'sxz'), 'exz': ('szz', 'sxz'), 'i ': ('tanPhi',)}
# better show rotation of particles
Gl1_Sphere.stripes = True
# open the plot on the screen
plot.plot()
O.saveTmp()
3d postprocessing¶
Following example is in file doc/sphinx/tutorial/05-3d-postprocessing.py. This example will run for 20000 iterations, saving *.png snapshots, then it will make a video 3d.mpeg out of those snapshots.
# demonstrate 3d postprocessing with yade
#
# 1. qt.SnapshotEngine saves images of the 3d view as it appears on the screen periodically
# makeVideo is then used to make real movie from those images
# 2. VTKRecorder saves data in files which can be opened with Paraview
# see the User's manual for an intro to Paraview
# generate loose packing
from yade import pack, qt
sp = pack.SpherePack()
sp.makeCloud((0, 0, 0), (2, 2, 2), rMean=.1, rRelFuzz=.6, periodic=True)
# add to scene, make it periodic
sp.toSimulation()
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
# interaction loop
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(damping=.4),
# save data for Paraview
VTKRecorder(fileName='3d-vtk-', recorders=['all'], iterPeriod=1000),
# save data from Yade's own 3d view
qt.SnapshotEngine(fileBase='3d-', iterPeriod=200, label='snapshot'),
# this engine will be called after 20000 steps, only once
PyRunner(command='finish()', iterPeriod=20000)
]
O.dt = .5 * PWaveTimeStep()
# prescribe constant-strain deformation of the cell
O.cell.velGrad = Matrix3(-.1, 0, 0, 0, -.1, 0, 0, 0, -.1)
# we must open the view explicitly (limitation of the qt.SnapshotEngine)
qt.View()
# this function is called when the simulation is finished
def finish():
# snapshot is label of qt.SnapshotEngine
# the 'snapshots' attribute contains list of all saved files
makeVideo(snapshot.snapshots, '3d.mpeg', fps=10, bps=10000)
O.pause()
# set parameters of the renderer, to show network chains rather than particles
# these settings are accessible from the Controller window, on the second tab ("Display") as well
rr = yade.qt.Renderer()
rr.shape = False
rr.intrPhys = True
Periodic triaxial test¶
Following example is in file doc/sphinx/tutorial/06-periodic-triaxial-test.py. A variant of this exemple includes capillary forces, see doc/sphinx/tutorial/06-periodic-triaxial-test-capillarity.py
# encoding: utf-8
# periodic triaxial test simulation
#
# The initial packing is either
#
# 1. random cloud with uniform distribution, or
# 2. cloud with specified granulometry (radii and percentages), or
# 3. cloud of clumps, i.e. rigid aggregates of several particles
#
# The triaxial consists of 2 stages:
#
# 1. isotropic compaction, until sigmaIso is reached in all directions;
# this stage is ended by calling compactionFinished()
# 2. constant-strain deformation along the z-axis, while maintaining
# constant stress (sigmaIso) laterally; this stage is ended by calling
# triaxFinished()
#
# Controlling of strain and stresses is performed via PeriTriaxController,
# of which parameters determine type of control and also stability
# condition (maxUnbalanced) so that the packing is considered stabilized
# and the stage is done.
#
sigmaIso = -1e5
#import matplotlib
#matplotlib.use('Agg')
# generate loose packing
from yade import pack, qt, plot
O.periodic = True
sp = pack.SpherePack()
if 0:
## uniform distribution
sp.makeCloud((0, 0, 0), (2, 2, 2), rMean=.1, rRelFuzz=.3, periodic=True)
else:
## create packing from clumps
# configuration of one clump
c1 = pack.SpherePack([((0, 0, 0), .03333), ((.03, 0, 0), .017), ((0, .03, 0), .017)])
# make cloud using the configuration c1 (there could c2, c3, ...; selection between them would be random)
sp.makeClumpCloud((0, 0, 0), (2, 2, 2), [c1], periodic=True, num=500)
# setup periodic boundary, insert the packing
sp.toSimulation()
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop([Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
PeriTriaxController(
label='triax',
# specify target values and whether they are strains or stresses
goal=(sigmaIso, sigmaIso, sigmaIso),
stressMask=7,
# type of servo-control
dynCell=True,
maxStrainRate=(10, 10, 10),
# wait until the unbalanced force goes below this value
maxUnbalanced=.1,
relStressTol=1e-3,
# call this function when goal is reached and the packing is stable
doneHook='compactionFinished()'
),
NewtonIntegrator(damping=.2),
PyRunner(command='addPlotData()', iterPeriod=100),
]
O.dt = .5 * PWaveTimeStep()
def addPlotData():
plot.addData(
unbalanced=unbalancedForce(),
i=O.iter,
sxx=triax.stress[0],
syy=triax.stress[1],
szz=triax.stress[2],
exx=triax.strain[0],
eyy=triax.strain[1],
ezz=triax.strain[2],
# save all available energy data
Etot=O.energy.total(),
**O.energy
)
# enable energy tracking in the code
O.trackEnergy = True
# define what to plot
plot.plots = {
'i': ('unbalanced',),
'i ': ('sxx', 'syy', 'szz'),
' i': ('exx', 'eyy', 'ezz'),
# energy plot
' i ': (O.energy.keys, None, 'Etot'),
}
# show the plot
plot.plot()
def compactionFinished():
# set the current cell configuration to be the reference one
O.cell.trsf = Matrix3.Identity
# change control type: keep constant confinement in x,y, 20% compression in z
triax.goal = (sigmaIso, sigmaIso, -.2)
triax.stressMask = 3
# allow faster deformation along x,y to better maintain stresses
triax.maxStrainRate = (1., 1., .1)
# next time, call triaxFinished instead of compactionFinished
triax.doneHook = 'triaxFinished()'
# do not wait for stabilization before calling triaxFinished
triax.maxUnbalanced = 10
def triaxFinished():
print('Finished')
O.pause()
Fluid injection¶
Following example is in file doc/sphinx/tutorial/07-fluid-injection.py. The video below results from post-processing with paraview
# This script simulates the injection of a fluid in a localized region below immersed particles
# The simulation is periodic along z-axis.
# at first execution, the simulation starts by depositing the particles in the container then saves the scene before proceeding to injection
# further execution will reload the deposited layer and start injection directly to gain time
# WARNING: changes in some input parameters like dimensions of the box or number of particles will not be reflected as long as the saved state is present on disk,
# remember to erase it to force a new generation, or set newSample=True below
from yade import pack, export
import yade.timing
from math import *
from pylab import rand
import os.path
import numpy as np
import matplotlib.pyplot as plt
O.periodic = True
# Dimensions of the box and of the injection zone
width = 0.8 #
height = 1
depth = 0.4
aperture = 0.05 * width
# number of spheres
numSpheres = 2000
# contact friction during deposition
compFricDegree = 10
# porosity of the initial cloud
porosity = 0.8
# Cundall's damping (zero recommanded)
damp = 0.0
# fluid viscosity
mu = 0.01
# flow rate at the inlet
injectedFlux = -0.001
# name of output folder
key = 'output0'
newSample = False #turn this true if you want to generate new sample by pluviation each time you run the script
# Deduced mean size for generating the cloud and consistency check
filename = "init" + key + str(numSpheres)
volume = width * height * depth
meanRad = pow(volume * (1 - porosity) / (pi * (4 / 3.) * numSpheres), 1 / 3.)
if (meanRad * 6 > depth):
print("INCOMPATIBLE SIZES. INCREASE DEPTH OR INCREASE NUM_SPHERES")
# if no deposited layer has been found (first execution), generate bodies
if not os.path.isfile(filename + ".yade") or newSample: #we create new sample if it does not exist...
O.cell.hSize = Matrix3(width, 0, 0, 0, 3. * height, 0, 0, 0, depth)
O.materials.append(FrictMat(young=400000.0, poisson=0.5, frictionAngle=compFricDegree / 180.0 * pi, density=1600, label='spheres'))
O.materials.append(FrictMat(young=400000.0, poisson=0.5, frictionAngle=radians(15), density=1000, label='walls'))
lowBox = box(center=(width / 2.0, height, width / 2.0), extents=(width * 1000.0, 0, width * 1000.0), fixed=True, wire=False)
O.bodies.append(lowBox)
topBox = box(center=(width / 2.0, 2 * height + 4 * meanRad, width / 2.0), extents=(width * 1000.0, 0, width * 1000.0), fixed=True, wire=False)
O.bodies.append(topBox)
sp = pack.SpherePack()
sp.makeCloud((0, height + 2 * meanRad, 0), (width, 2 * height + 2 * meanRad, depth), -1, .002, numSpheres, periodic=True, porosity=porosity, seed=2)
O.bodies.append([sphere(s[0], s[1], color=(0.6 + 0.15 * rand(), 0.5 + 0.15 * rand(), 0.15 + 0.15 * rand()), material='spheres') for s in sp])
else: #... else we re-use the previous one
O.load(filename + ".yade")
# look better
Gl1_Sphere.stripes = True
# define the fluid solver with appropriate parameter (see PeriodicFlowEngine's documentation)
flow = PeriodicFlowEngine(
dead=1,
meshUpdateInterval=40,
defTolerance=-1,
permeabilityFactor=1.0,
useSolver=3,
duplicateThreshold=depth,
wallIds=[-1, -1, 0, 1, -1, -1],
bndCondIsPressure=[0, 0, 0, 1, 0, 0],
viscosity=mu,
label="flow"
)
newton = NewtonIntegrator(damping=damp, gravity=(0, -9.81, 0))
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Box_Aabb()], allowBiggerThanPeriod=1),
InteractionLoop([Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
GlobalStiffnessTimeStepper(active=1, timeStepUpdateInterval=1000, timestepSafetyCoefficient=0.3, defaultDt=utils.PWaveTimeStep(), label='timestepper'),
flow, # <====== the solver is inserted here, for the moment it is "dead" (doing nothing)
newton,
# some recorders for post-processing
PyRunner(command="flow.saveVtk(key)", iterPeriod=25, dead=1, label="vtkE"),
VTKRecorder(recorders=["spheres", "velocity", "stress"], iterPeriod=25, dead=1, fileName=key + '/', label="vtkR")
]
########## if this is fresh execution, get static equilibrium and save the result for later use #############
if not os.path.isfile(filename + ".yade") or newSample:
O.run(1000, 1)
while unbalancedForce() > 0.01:
O.run(100, 1)
# turn the recorders and the solver on
vtkE.dead = vtkR.dead = flow.dead = 0
#add a recorder and define what to plot
O.engines = O.engines + [PyRunner(command=('myAddPlotData()'), iterPeriod=50, label="recorder")]
O.save(filename + ".yade")
O.saveTmp()
########## define what to plot ##############
from yade import plot
import pylab
# First, find particle at the top of the sample (for evaluating initial height of the layer)
maxY = 0
for s in O.bodies:
if isinstance(s.shape, Sphere):
pos = s.state.pos
if pos[1] > maxY:
maxY = pos[1]
def myAddPlotData():
index = flow.getCell(0.5 * width, height, 0.5 * depth)
if index > 0:
########## find particle at the top of the sample ##############
simpleH = 0
for s in O.bodies:
if isinstance(s.shape, Sphere):
pos = s.state.pos
if pos[1] > simpleH:
simpleH = pos[1]
########## function to compute hf ##############
cavityh = height
for s in O.bodies:
v = s.state.vel
magvel = pow((v[0] * v[0] + v[1] * v[1] + v[2] * v[2]), 0.5)
if magvel > 0.14:
pos = s.state.pos
if pos[1] > cavityh:
cavityh = pos[1]
#########################
plot.addData(
t=O.time,
i=O.iter,
p=flow.getPorePressure((0.5 * width, height, 0.5 * depth)),
q=-injectedFlux,
Ho=maxY - 1,
hf=cavityh - 1,
H=simpleH - 1
)
H = simpleH - 1
hf = cavityh - 1
################ impose the costant flux ###############
# In this function we find the elements of the mesh which have a face in the injection region, to distribute the inlet flux
# it is inserted in the solver below
def imposeFlux(valF):
found = 0
listF1 = []
for k in range(flow.nCells()):
if 0 in flow.getVertices(k) and flow.getCellCenter(k)[0] > ((width - aperture) / 2.) and flow.getCellCenter(k)[0] < ((width + aperture) / 2.):
listF1.append(k)
flow.clearImposedFlux()
if len(listF1) == 0:
flow.imposeFlux((0.5 * width, height, 0.5 * depth), valF)
else:
for k in listF1:
flow.imposeFlux(flow.getCellBarycenter(k), valF / float(len(listF1)))
injectedFlux = -1 #this one is large and it will fluidize violently, see below for smoother evolutions
# very important: the flux needs to be imposed with this "hook", which plugs our custom function in the right place in the solving sequence
flow.blockHook = "imposeFlux(injectedFlux)"
############### EXERCISE ###################
# 1- use trial and error in the shell (see the note below) to find approximately an upper-bound of "injectedFlux", below which there is only limited movements of the particles
# 2- implement in this script a progressive increase of the flux, starting from a small fraction of max value above, and exceeding it by a factor 2 at the end
# 3- choose an approriate plot to show that the pressure response is initially linear, then sublinear. Is the upper-bound from question 1 also an upper-bound for the linear response?
# 4- compare with fig. 10 from https://doi.org/10.1103/PhysRevE.94.052905 (also available at https://arxiv.org/pdf/1703.02319)
# 5- use paraview to generate a video similar to https://www.youtube.com/watch?v=gH585XaQEcY (it can be with a constant flux)
#NOTE:
# to change the injected flux interactively, change it in the global scope:
# globals()["injectedFlux"]=-0.03
# else we have two variables with the same name in different scopes and
