New question #680746 on Yade:
https://answers.launchpad.net/yade/+question/680746
Hi everyone,
I am new to Yade and python and I am trying to simulate a pile and a shallow as
a clump. One friend helps me to create a script to generate a sample of my soil
and my foundation. I'm studying geotechinical problems. But i don't know how
to obtain the total Force acting in the clump, because of a velocity imposed.
Can someone tell me how to do this? This is my script:
# -*- coding: utf-8 -*-
from yade import pack,utils,plot,ymport
import math
################
### 1.INPUTS ###
################
# 1.1 IMPORT PARAMETERS FROM TABLE
nRead=readParamsFromTable(
DensidadeP = 2.70,
DensidadeS = 3.50,
YoungP = 1200e6,
PoissonP = 0.20,
AnguloAtritoP = 39.5, #variar angulo de atrito para calibrar com a
curva do professor Nakai
Porosidade = 0.20,
MatType = 'frict',
TestDim = '2D',
rate = 0.1,
PSD = [[0.0032,0.006],[0.6,1.]],
TensaoIsotropica = 1.9e4,
verlet = 0.035,
unknownOk=True
)
from yade.params import table
# 1.2 PARTICLE/MEDIUM PARAMETERS
mediumPorosity = table.Porosidade
particleDensity = table.DensidadeP #kg/m³
particleYoung = table.YoungP #Pa
particlePoisson = table.PoissonP #Adimensional
particleFricAng = radians(table.AnguloAtritoP)
# 1.3 ELEMENTS DIMENSIONS
# 1.3.1 MEDIUM PARTICLE SIZE DISTRIBUTION
PSD = table.PSD
# 1.3.2 SPHERES RADIUS
meanElementRadus = ((PSD[0][0]+PSD[0][-1])/2)/2
maxElementRadus=PSD[0][-1]/2
minElementRadus=PSD[0][0]/2
# 1.4 MATERIAL TYPE
materialType = 'frict' # 'frict', 'cohfrict', 'concrete'
# 1.5 DIMENSIONS (2D OR 3D)
testDimensions = table.TestDim
# 1.6 SHALLOW AND PILE (2D)
AlturaBloco=0.03 #Altura do
bloco de fundação
BordasBloco=0.055 #Exedente dos bordos do
bloco com relação à(s) estaca(s)
ProfEstaca=0.24 #Profundidade da estaca
CompPonta=0 #Comprimento da ponta
da estaca (estacas com ponta ou tubulões)
DiametroEstaca=0.01 #Diâmetro da estaca
Afastamentos=5*DiametroEstaca #Região do solo nao
influenciada pela estaca
XCentroE=0. #Coordenada 'x' do
centro da estaca
YCentroE1=0.
YCentroE2=0.04
YCentroE3=-0.04 #Coordenada 'y' do
centro da estaca
ZCentroE= (0.5-ProfEstaca/2) #Coordenada 'z' do centro da
estaca
CentroEstaca1=(XCentroE,YCentroE1,ZCentroE) #Vetor de coordenadas
do centro da estaca
TopoEstaca1=(XCentroE,YCentroE1,ZCentroE-ProfEstaca/2) #Vetor de coordenadas
do topo da estaca
BaseEstaca1=(XCentroE,YCentroE1,ZCentroE+ProfEstaca/2) #Vetor de coordenadas
da ponta da estaca
CentroEstaca2=(XCentroE,YCentroE2,ZCentroE) #Vetor de coordenadas
do centro da estaca
TopoEstaca2=(XCentroE,YCentroE2,ZCentroE-ProfEstaca/2) #Vetor de coordenadas
do topo da estaca
BaseEstaca2=(XCentroE,YCentroE2,ZCentroE+ProfEstaca/2) #Vetor de coordenadas
da ponta da estaca
TopoEstaca3=(XCentroE,YCentroE3,ZCentroE-ProfEstaca/2) #Vetor de coordenadas
do topo da estaca
BaseEstaca3=(XCentroE,YCentroE3,ZCentroE+ProfEstaca/2)
XCentroB=XCentroE #Coordenada 'x' do
centro do bloco
YCentroB=YCentroE1 #Coordenada 'y' do
centro do bloco
ZCentroB=0.5+AlturaBloco/2 #Coordenada 'z'
do centro do bloco
CentroBloco=(XCentroB,YCentroB,ZCentroB) #Vetor de coordenadas
do centro do bloco
Vertice1Bloco=(-1.01*1*0.003,YCentroB-DiametroEstaca-BordasBloco,ZCentroB-AlturaBloco/2)
Vertice2Bloco=(1.01*1*0.003,YCentroB+DiametroEstaca+BordasBloco,ZCentroB+AlturaBloco/2)
# 1.7 IMPORT/EXPORT NAMES
ImportName = 'Packing' + testDimensions + str(materialType)+ '-' +
str(table.Porosidade)
ExportName = 'Triaxial' + testDimensions + str(materialType) + '-Phi' +
str(table.AnguloAtritoP) + '-Y' + str(particleYoung/1.e6)+'MPa'+ '-Poisson' +
str(particlePoisson) + '-' + str(mediumPorosity) + ' - Iso' +
str(table.TensaoIsotropica)
# 1.8. USER DEFINED FUNCTIONS
# 1.8.1. POROSITY
def Porosity(Region):
soma=0.
for b in O.bodies:
if isinstance(b.shape,Sphere):
soma+=pi*(b.shape.radius**2)
return (Region-soma)/Region
################## ### 2.MATERIAL ### ##################
# 2.1 MEDIUM MATERIAL
if materialType == 'frict':
spheresMaterial=FrictMat(density=particleDensity,
frictionAngle=particleFricAng, young=particleYoung, poisson=particlePoisson)
O.materials.append(spheresMaterial)
# 2.2 BOUNDARY MATERIAL
boundaryMaterial=FrictMat(density=particleDensity,
frictionAngle=particleFricAng, young=10.*particleYoung,
poisson=2*particlePoisson)
O.materials.append(boundaryMaterial)
#2.3 SHALLOW AND PILE MATERIAL
AngAtritoA=atan(0)
DensidadeA=27000.
YoungA=70.e9
PoissonA=.33
MatConcreto=O.materials.append(FrictMat(density=DensidadeA,young=YoungA,poisson=PoissonA,frictionAngle=AngAtritoA,label="aluminio"))
################## ### 3.GEOMETRY ### ##################
# 3.1 PACKING DIMENSIONS
width = 1
height = 0.5
dx = 1.01*2*maxElementRadus
dy = width
dz = height
CorteTransversal = dy
baseArea = dy
lateralAreaX = dy
lateralAreaY = 1.
# 3.2 BOUNDARIES
corner1,corner2 = (-.5*dx,-.5*dy,0.),(.5*dx,.5*dy,dz)
walls = aabbWalls([corner1, corner2], thickness=0., oversizeFactor=2,
material=boundaryMaterial)
wallIds = O.bodies.append(walls)
base=O.bodies[4]
topo=O.bodies[5]
#3.3 SHALLOW AND PILE
GeomEst1=pack.inCylinder(BaseEstaca1,TopoEstaca1,DiametroEstaca/2)
GeomEst2=pack.inCylinder(BaseEstaca2,TopoEstaca2,DiametroEstaca/2)
GeomEst3=pack.inCylinder(BaseEstaca3,TopoEstaca3,DiametroEstaca/2)
GeomBloco=pack.inAlignedBox(Vertice1Bloco,Vertice2Bloco)
GeomBloco=pack.regularOrtho(GeomBloco,radius=0.002,gap=0.,color=(0,0,1),material="aluminio")
GeomEst1=pack.regularOrtho(GeomEst1,radius=0.002,gap=0.,color=(0,0,1),material="aluminio")
GeomEst2=pack.regularOrtho(GeomEst2,radius=0.002,gap=0.,color=(0,0,1),material="aluminio")
GeomEst3=pack.regularOrtho(GeomEst3,radius=0.002,gap=0.,color=(0,0,1),material="aluminio")
O.bodies.appendClumped(GeomBloco)
O.bodies.appendClumped(GeomEst1)
O.bodies.appendClumped(GeomEst2)
O.bodies.appendClumped(GeomEst3)
# 3.4 IMPORT SPHERE PACKING
O.bodies.append(ymport.textExt(fileName = ImportName, format='x_y_z_r_matId'))
# 3.5. ASSIGN MATERIAL PROPERTIES TO THE SPHERES
for b in O.bodies:
if isinstance(b.shape,Sphere):
b.material=spheresMaterial
# 3.6. DEFINE BLOCKED DOFs
for b in O.bodies:
if isinstance(b.shape,Sphere):
b.state.blockedDOFs='xYZ'
elif b.isClump:
b.state.blockedDOFs='xXYZ'
# 4.1. TEST PROCEDURE
def Analise():
O.forces.f(topo.id)[2]
def Porosity(Region):
if testDimensions == '3D':
return porosity(Region)
elif testDimensions == '2D':
soma=0.
for b in O.bodies:
if isinstance(b.shape,Sphere):
soma+=pi*(b.shape.radius**2)
return (Region-soma)/Region
if testDimensions == '2D':
acelerador = 2
elif testDimensions == '3D':
acelerador = 5
topo.state.pos[2] = topo.state.pos[2]-0.24
def PreencherTriaxial():
if Porosity(CorteTransversal*topo.state.pos[2])>0.13:
for b in O.bodies:
if b.isClump:
b.state.vel = (0,0,-10*acelerador)
newton.damping = 0.7
else:
print 'fim'
#def force():
# global force
# if O.iter>0.:
# i=6.
# forca=0.
# while i>5 and i<352:
# forca= forca +O.forces.f(i)[2]
# i=i+1
# return list(force)
def myAddData():
b=O.bodies[351]
plot.addData(z1=0.53-b.state.pos[2], f1=O.forces.f(6)[2], i=O.iter,
t=O.time)
plot.plots={'t':('f1','z1'),'f1':'z1'}
plot.plot()
O.run()
utils.waitIfBatch()
############################### ### 5. RECORD AND PLOT DATA ###
###############################
#def Displacement():
# global Displacement
# for b in b.isClump:
# print b.state.pos[2]
############################## ### 6.SIMULATION PROCEDURE ###
##############################
enlargeFactor = 1.5
O.engines=[
ForceResetter(),
InsertionSortCollider([
Bo1_Sphere_Aabb(),
Bo1_Box_Aabb(),
Bo1_Wall_Aabb(),
Bo1_Facet_Aabb()
],verletDist=table.verlet*2.*meanElementRadus),
InteractionLoop(
[Ig2_Wall_Sphere_ScGeom(),
Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor),
Ig2_Box_Sphere_ScGeom(),
Ig2_Facet_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_MindlinPhys(en=0.8, es=0.8)],
[Law2_ScGeom_MindlinPhys_Mindlin()],
),
GlobalStiffnessTimeStepper(
active=True,
defaultDt=SpherePWaveTimeStep(radius=meanElementRadus,
density=O.materials[0].density, young=O.materials[0].young),
timeStepUpdateInterval=1000,timestepSafetyCoefficient=0.5
),
# PyRunner(command='force()',iterPeriod=10,label='PrepararAmostra'),
PyRunner(command='PreencherTriaxial()',iterPeriod=10,label='PrepararAmostra'),
PyRunner(command='myAddData()',iterPeriod=10,label='PrepararAmostra'),
NewtonIntegrator(damping=0.7, gravity=[0.,0.,-9.81], label='newton'),
]
########################## ### 7. RUN SIMULATION ###
##########################
O.run()
utils.waitIfBatch()
--
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
