how to set flexible boundaries?

Asked by Xifan Li

I want to set the flexible boundary in Yade.

From the answer https://answers.launchpad.net/yade/+question/254339
I found that there are several options to achieve this. Such as the bonded particles, and the flexible membrane by woo.

But I want to try another way:
[1] finish the same task in Fig.9 as this paper: https://ascelibrary.org/doi/full/10.1061/(ASCE)GM.1943-5622.0001001
[2] In the above paper, the author set the boundary as the bonded particles, But I want to import velocity (stress) on each particle.
[3]Maybe we can try to set the moving rigid boundary so that So we can look at it approximately as a flexible boundary?

I do not know which way is the best and achievable.

Any reply would be a grateful appreciate.

My Yade code is as followed, which was modified from the Oedemetric Test from the Examples, and I do not know how to modify it.
#################################################
# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(rMean=.03,rRelFuzz=0,maxLoad=1.5e4,minLoad=1e1)
# 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((1,0,1),(1,.03,1),wallMask=31))
#O.bodies.append(geom.facetBox((0,0,0),(.5,.03,.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(2,0.03,2),rMean=rMean,rRelFuzz=rRelFuzz)
#sp.makeCloud((0,0,0),(-0.5,0.03,-0.25),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=1,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, -.05)
 # start plotting the data now, it was not interesting before
 O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=50)]
 # 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()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Xifan Li
Solved:
Last query:
Last reply:
Revision history for this message
Xifan Li (shelvan) said (last edit ):
#1

I have tried to set the confining pressure into this model. But it didn't work. An error occurred while running.
I added the TriaxialStressController into the previous code.

The following is my modified code:

from __future__ import print_function
sigmaIso = -1e5

# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(rMean=.03,rRelFuzz=0,maxLoad=1.5e4,minLoad=1e1)
# 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, qt, plot

O.periodic = True

O.bodies.append(geom.facetBox((2,0,2),(2,.03,2),wallMask=31))
#O.bodies.append(geom.facetBox((0,0,0),(.5,.03,.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((1,0,0),(3,0.03,3.5),rMean=rMean,rRelFuzz=rRelFuzz, periodic = True)
#sp.makeCloud((0,0,0),(-0.5,0.03,-0.25),rMean=rMean,rRelFuzz=rRelFuzz)
sp.toSimulation()

triax=TriaxialStressController(
 ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.

 thickness = 0,
 stressMask = 2,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

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()]),
  PeriTriaxController(
                label='triax',
                # specify target values and whether they are strains or stresses
                goal=(sigmaIso, sigmaIso, sigmaIso),
                stressMask=2,
                # type of servo-control
                dynCell=True,
                maxStrainRate=(10, 10, 10),
                # wait until the unbalanced force goes below this value
                maxUnbalanced=10,
                relStressTol=1e-3,

        ),

 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=1,label='checker'),
]
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

triax.goal1=triax.goal2=triax.goal3=-10000

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, -.05)
 # start plotting the data now, it was not interesting before
 O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=50)]
 # 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 *= -.5
  # 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()

def triaxFinished():
 print('Finished')
 O.pause()

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()
#######################################################

Could you please help me to check about my code and tell me which part goes wrong?

Revision history for this message
Jan Stránský (honzik) said :
#2

Hello,

> I found that there are several options to achieve this

yes, many options.
Which one is suitable depends on many factors.

> [2] In the above paper, the author set the boundary as the bonded particles, But I want to import velocity (stress) on each particle.

you can set velocity of any particle [2] or force to any particle (simulating stress) [3].

> [3]Maybe we can try to set the moving rigid boundary so that So we can look at it approximately as a flexible boundary?

Maybe.
This you should determine, what is suitable for your task.

> I do not know which way is the best and achievable.

nobody knows.. it strongly depends on definition of "best" and "achievable".

> But it didn't work. An error occurred while running.

please read [1] and be more specific

> Could you please help me to check about my code and tell me which part goes wrong?

Please try to make the code a MWE first [1]

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.State.vel
[3] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ForceContainer.addF

Revision history for this message
Xifan Li (shelvan) said :
#3