Question #679808 on Yade changed:
https://answers.launchpad.net/yade/+question/679808
Description changed to:
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 -*-
from yade import pack,utils,plot,export
import numpy as np
import operator,fractions
# TODAS AS UNIDADES ESTÃO NO S.I.
###############################
### 1.PARÂMETROS DE ENTRADA ###
###############################
# 1.1 IMPORT PARAMETERS FROM TABLE
nRead=readParamsFromTable(
DensidadeP = 2.7,
PoissonP = 0.2,
AnguloAtritoP = 0,
CoesaoP = 0,
PorosidadeS = 0.20,
MatType = 'frict',
TestType = '2D',
PSD = [[0.0032,0.006],[0.6,1.]],
verlet = 0.35,
YoungP = 7.e4,
unknownOk=True
)
from yade.params import table
# 1.2 PARAMETROS DO MEIO (MACROSCÓPICOS)
DensidadeS = table.DensidadeP #tf/m³
YoungP = table.YoungP #kN/m²
PoissonP = table.PoissonP #Adimensional
AnguloAtritoP = radians(table.AnguloAtritoP)
CoesaoP = table.CoesaoP
PorosidadeS = table.PorosidadeS
# 1.3 DIMENSÕES DOS ELEMENTOS
# 1.3.1 CURVA GRANULOMETRICA DO SOLO - PSD=[[diâmetros],[Quantidade
passante Acumulada]]
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,
defaultDt=SpherePWaveTimeStep(radius=rMeanElemento,
density=O.materials[0].density, young=O.materials[0].young),
timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8
),
NewtonIntegrator(damping=0.3, gravity=[0.,0.,-9.81], label='newton'),
PyRunner(command='PreencherTriaxial()',iterPeriod=10,label='PrepararAmostra'),
# PyRunner(command='addPlotData()',iterPeriod=10),
]
##################
### 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):
soma+=pi*(b.shape.radius**2)
return (Region-soma)/Region
if TestType == '2D':
acelerador = 1
elif TestType == '3D':
acelerador = 5
Compressao=False
Descompressao=False
AlturaPreenchida=0
def PreencherTriaxial():
global Compressao
global Descompressao
global topo
global AlturaPreenchida
if Compressao:
if Porosity(CorteTransversal*topo.state.pos[2])<PorosidadeS:
topo.state.vel = (0.,0.,1.*acelerador)
Descompressao = True
Compressao = False
newton.damping = 0.7
else:
topo.state.vel = (0,0,-10*acelerador)
elif Descompressao:
if abs(O.forces.f(topo.id)[2])==0:
if topo.state.vel==(0,0,1*acelerador):
print
Porosity(CorteTransversal*topo.state.pos[2])
AlturaPreenchida=topo.state.pos[2]
print AlturaPreenchida/height
if topo.state.pos==(0,0,dz+rMaxElemento):
Material.density=DensidadeS*100/(1-Porosity(CorteTransversal*AlturaPreenchida))
if AlturaPreenchida>height-rMeanElemento:
topo.state.vel=(0,0,0)
calm()
Descompressao=False
Material.frictionAngle=radians(35)
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,30*acelerador)
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 "Amostra Preparada"
print "Densidade"
print Material.density
print "Porosidade"
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.
--
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 : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
