Question #699605 on Yade changed:
https://answers.launchpad.net/yade/+question/699605
Status: Needs information => Open
Fedor gave more information on the question:
Thank you for notes. My apologies
I want to use orthogonal packing and investigate how to change porosity vs
isotropic pressure with diffrent values of young modulus.
I used the existing script and added young modulus, but when I vary module the
porosity changes inadequate. There is the script.
from __future__ import print_function
O.periodic = True
O.cell.hSize = Matrix3(2, 0, 0, 0, 2, 0, 0, 0, 2)
from yade import pack, plot
sigmaIso = -7e5
young=5e8
poisson=0.2
frictionAngle=0
# 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
# in this case, add dense packing
pp=O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=frictionAngle))
O.bodies.append(pack.regularOrtho(pack.inAlignedBox((0, 0, 0), (2, 2, 2)),
radius=.09, gap=0, color=(0, 0, 1), material=pp))
#O.bodies.append(pack.regularHexa(pack.inAlignedBox((0, 0, 0), (2, 2,
2)), radius=.09, 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),
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=5,
relStressTol=1e-3,
# call this function when goal is reached and the packing is
stable
doneHook='compactionFinished()'
),
]
# 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(0, 0, 0, 0, -.1, 0, 0, 0, 0)
# 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)
porrr=utils.voxelPorosity(200,(0.75,0.75,0.75),(1.25,1.25,1.25))
#porrr=utils.aabbExtrema(0, False)
print('porosity', porrr)
# called from the 'checker' engine periodically, during the shear phase
# 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,
porrr=utils.voxelPorosity(200,(0.75,0.75,0.75),(1.25,1.25,1.25)), stress =
getStress().trace() / 3.)
# 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 ': ('stress',), 'stress':('porrr')}
# better show rotation of particles
Gl1_Sphere.stripes = True
# open the plot on the screen
plot.plot()
O.saveTmp()
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, sigmaIso)
triax.stressMask = 7
# 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()
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
