yade-users team mailing list archive

[Fedor <question699605@xxxxxxxxxxxxxxxxxxxxx>]

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.