# yade-users team mailing list archive

## [Question #679808]: rigid boundaries

```New question #679808 on Yade:

Hi everyone:

I am new to Yade and python and I am trying to simulate a pile as a rigid boundaries. One friend helps me to create a script to generate a sample of my soil. I'm studying geotechinical problems.  But i don't know how to create a new rigid boundary to simulate a pile. Can someone tell me how to do this? This is my script:

# -*- coding: utf-8 -*-

#In order to match simulations to experiments, a procedure to calibrate these parameters must be given.
#Since elastic properties are independent of inelastic ones, they are calibrated first.

#Order of calibration:
#	1. elastic properties, which depend on only geometry and elastic parameters
#		E - Youngs modulus,
#		ν - Poisson’s ratio
#	2. inelastic properties, depending (in addition) on damage and plasticity parameters:
#		ft - tensile strength
#		fc - compressive strength
#		Gf - fracture energy
#	3. confinement properties; they appear only in high confinement situations and can be calibrated without having substantial impact on
#	   already calibrated inelastic properties. We do not describe them quantitatively; fitting simulation and experimental curves is used
#	4. rate-dependence properties; they appear only in high-rate situations, therefore are again calibrated after inelastic properties
#	   independently. As in the previous case, a simple fitting approach is used here.

### PRocedimento:
# 0. Encher a caixa - Bernardo
# 1. Calibrar densidade - Bernardo
# 2. Young
#	- Bruna:
#		-Pesquisar obtenção do parametro pelo ensaio triaxial. Ver se a tensão desviadora influencia o processo.
#	-Bernardo:
# 		-Implementar o cálculo
# 3. Poisson
#	-Como implementar isso? Como obter ele do ensaio triaxial??
# 4. Angulo de Atrito
#	- Ensaio de cisalhamento direto ou triaxial??
# 5. Coesão

import numpy as np
import operator,fractions

# TODAS AS UNIDADES ESTÃO NO S.I.

###############################
###############################

# 1.1 IMPORT PARAMETERS FROM TABLE
PoissonP = 0.2,
AnguloAtritoP = 0,
CoesaoP = 0,
MatType = 'frict',
TestType = '2D',
PSD = [[0.0032,0.006],[0.6,1.]],
verlet = 0.35,
YoungP = 7.e4,
unknownOk=True
)

# 1.2 PARAMETROS DO MEIO (MACROSCÓPICOS)
YoungP = table.YoungP	#kN/m²
CoesaoP = table.CoesaoP

# 1.3 DIMENSÕES DOS ELEMENTOS

PSD = table.PSD

# 1.3.2 DADOS MÉDIOS DOS ELEMENTOS
rMeanElemento = ((PSD[0][0]+PSD[0][-1]))/2
rRelFuzzElemento = PSD[0][-1]/(2*rMeanElemento)-1
rMaxElemento = PSD[0][-1]/2
rMinElemento = PSD[0][0]/2

# 1.4 TIPO DE MATERIAL
MatType = table.MatType 	# 'frict', 'cohfrict', 'concrete'

# 1.5 TIPO DE AMOSTRA
TestType = table.TestType

##################
### 2.MATERIAL ###
##################

# 2.1 Material Meio
if MatType == 'frict':
Material=FrictMat(density=20400, frictionAngle=0)
elif MatType == 'cohfrict':
Material=CohFrictMat(density=2600)
elif MatType == 'concrete':
#	Material=CpmMat(density=2600)
Material=FrictMat(density=2600, frictionAngle=0)
O.materials.append(Material)

# 2.2 Material Contorno
MatContorno = FrictMat(frictionAngle=0)
O.materials.append(MatContorno)

###################
### 3.GEOMETRIA ###
###################

# 3.1 Triaxial

# 3.1.1 Dimensões
width = 1
height = 0.5

if TestType == '2D':
dx = 1.01*2*rMaxElemento
dy = width
dz = height
CorteTransversal = dy
elif TestType == '3D':
dx = width
dy = width
dz = height
CorteTransversal = dy*dx

# 3.1.2 Faces
mn,mx=(-.5*dx,-.5*dy,0),(.5*dx,.5*dy,dz)
walls=aabbWalls([mn,mx],thickness=0,material=MatContorno)
wallIds=O.bodies.append(walls)

base=O.bodies[4]
topo=O.bodies[5]

##################
### 4.CÁLCULOS ###
##################

# 4.1. InteractionLoop
enlargeFactor=1.5

LoopFrict=InteractionLoop(
[
Ig2_Wall_Sphere_ScGeom(),
Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor),
Ig2_Box_Sphere_ScGeom(),
Ig2_Facet_Sphere_ScGeom(),
],
[
Ip2_FrictMat_FrictMat_FrictPhys(),
],
[
Law2_ScGeom_FrictPhys_CundallStrack(),
],
)

LoopCohFrict=InteractionLoop(
[
Ig2_Wall_Sphere_ScGeom(),
Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=enlargeFactor),
Ig2_Box_Sphere_ScGeom6D(),
Ig2_Facet_Sphere_ScGeom6D(),
],
[
Ip2_FrictMat_FrictMat_FrictPhys(),
Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True,label='ipf')
],
[
Law2_ScGeom_FrictPhys_CundallStrack(),
Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True)
],
)
LoopConcrete=InteractionLoop(
[
Ig2_Wall_Sphere_ScGeom(),
Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor),
Ig2_Box_Sphere_ScGeom(),
Ig2_Facet_Sphere_ScGeom(),
],
[
Ip2_FrictMat_FrictMat_FrictPhys(),
Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=10),
Ip2_FrictMat_CpmMat_FrictPhys(),
],
[
Law2_ScGeom_FrictPhys_CundallStrack(),
Law2_ScGeom_CpmPhys_Cpm(),
],
)

if MatType == 'frict':
IntLoop = LoopFrict
elif MatType == 'cohfrict':
IntLoop = LoopCohFrict
elif MatType == 'concrete':
IntLoop = LoopConcrete

# 4.2 Factory Parameters

MassFlowRate = 5000/(CorteTransversal*dz)
VMAX = 5
VMIN = 5
FCenter = (0.,0.,height/2)
FExtents2D = (0,width/2.,height/2)
FExtents3D = (width/2.,width/2.,height/2)

factoy2D = BoxFactory(
center = FCenter, extents = FExtents2D,
PSDsizes = PSD[0], PSDcum = PSD[1], PSDcalculateMass = True,
vMax = VMAX, vMin = VMIN, vAngle = 0., normal = (0.,0.,-1.),
maxParticles = -1,
massFlowRate = MassFlowRate,
#		exactDiam = False,
label = 'factory',
blockedDOFs = 'xYZ',
materialId = 0
)

factoy3D = BoxFactory(
center = FCenter, extents = FExtents3D,
PSDsizes = PSD[0], PSDcum = PSD[1], PSDcalculateMass = True,
vMax = 5*VMAX, vMin = 5*VMIN, vAngle = 0., normal = (0.,0.,-1.),
maxParticles = -1,
massFlowRate = MassFlowRate,
#		exactDiam = False,
label = 'factory',
materialId = 0
)

if TestType == '2D':
factory = factoy2D
elif TestType == '3D':
factory = factoy3D

#4.3 Engines
O.engines=[
ForceResetter(),
InsertionSortCollider([
Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor),
Bo1_Box_Aabb(),
Bo1_Wall_Aabb(),
Bo1_Facet_Aabb()
], verletDist=0.07*2*rMeanElemento),
IntLoop,
factory,
DomainLimiter(lo=(-dx/2,-dy/2.,0.),hi=(dx/2,dy/2.,dz+5*rMeanElemento),iterPeriod=100),
GlobalStiffnessTimeStepper(
active=True,
timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8
),
NewtonIntegrator(damping=0.3, gravity=[0.,0.,-9.81], label='newton'),
PyRunner(command='PreencherTriaxial()',iterPeriod=10,label='PrepararAmostra'),
]

##################
### 5.PyRunner ###
##################

def Porosity(Region):
if TestType == '3D':
return porosity(Region)
elif TestType == '2D':
soma=0.
for b in O.bodies:
if isinstance(b.shape,Sphere):
return (Region-soma)/Region

if TestType == '2D':
elif TestType == '3D':

Compressao=False
Descompressao=False
AlturaPreenchida=0

def PreencherTriaxial():
global Compressao
global Descompressao
global topo
global AlturaPreenchida
if Compressao:
Descompressao = True
Compressao = False
newton.damping = 0.7
else:
elif Descompressao:
if abs(O.forces.f(topo.id)[2])==0:
print Porosity(CorteTransversal*topo.state.pos[2])
AlturaPreenchida=topo.state.pos[2]
print AlturaPreenchida/height
if topo.state.pos==(0,0,dz+rMaxElemento):
if AlturaPreenchida>height-rMeanElemento:
topo.state.vel=(0,0,0)
calm()
Descompressao=False
newton.damping=0.4
PrepararAmostra.command='EstabilizarTriaxial()'
PrepararAmostra.iterPeriod=1
else:
topo.state.vel = (0,0,0)
Descompressao = False
factory.massFlowRate = 220/(CorteTransversal*(height-AlturaPreenchida))
#					factory.center[2] = (AlturaPreenchida+height)/2
#					factory.extents[2] = (height-AlturaPreenchida)/2
newton.damping = 0.1
#					if AlturaPreenchida > 0.95*height:
#						newton.gravity = [0.,0.,-9.81]
else:
topo.state.vel=(0,0,0)
topo.state.pos=(0,0,dz+rMaxElemento)
elif factory.massFlowRate==0:
Compressao=True

def EstabilizarTriaxial():
if topo.state.pos[2]>0.98*height:
topo.state.vel=(0,0,-25)
else:
topo.state.vel=(0,0,0)
if utils.unbalancedForce()<0.005:      #Checar isso para concreto
PrepararAmostra.command='FinalizarPreparo()'

def FinalizarPreparo():
print Material.density
print Porosity(CorteTransversal*AlturaPreenchida)
print "Tempo"
print O.realtime
export.textExt('Packing' + TestType + MatType + '-' + str(PorosidadeS),'x_y_z_r_matId')
O.pause()

O.run()
utils.waitIfBatch()

Thank you so much and wish have a good day!

Regards,
Jessica.

--