# How to plot the total Force on a Clump

Asked by Jessica Soares da Rocha on 2019-05-10

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 -*-

import math

################
### 1.INPUTS ###
################

# 1.1 IMPORT PARAMETERS FROM TABLE
YoungP = 1200e6,
PoissonP = 0.20,
AnguloAtritoP = 39.5, #variar angulo de atrito para calibrar com a curva do professor Nakai
MatType = 'frict',
TestDim = '2D',
rate = 0.1,
PSD = [[0.0032,0.006],[0.6,1.]],
TensaoIsotropica = 1.9e4,
verlet = 0.035,
unknownOk=True
)

# 1.2 PARTICLE/MEDIUM PARAMETERS
particleYoung = table.YoungP #Pa

# 1.3 ELEMENTS DIMENSIONS

# 1.3.1 MEDIUM PARTICLE SIZE DISTRIBUTION
PSD = table.PSD

# 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)
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):
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)
YoungA=70.e9
PoissonA=.33

################## ### 3.GEOMETRY ### ##################

# 3.1 PACKING DIMENSIONS

width = 1
height = 0.5
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)

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):
return (Region-soma)/Region

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

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:
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)

b=O.bodies[351]

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()
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,
timeStepUpdateInterval=1000,timestepSafetyCoefficient=0.5
),
# PyRunner(command='force()',iterPeriod=10,label='PrepararAmostra'),
PyRunner(command='PreencherTriaxial()',iterPeriod=10,label='PrepararAmostra'),
NewtonIntegrator(damping=0.7, gravity=[0.,0.,-9.81], label='newton'),
]

########################## ### 7. RUN SIMULATION ### ##########################

O.run()
utils.waitIfBatch()

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
2019-05-13
Last query:
2019-05-13
2019-05-13
 Bruno Chareyre (bruno-chareyre) said on 2019-05-12: #1

Hi,
I don't think clump or imposed velocity make a difference:
O.forces.f[id] should be the force.

If you need help for a specific script please make it minimal. Nobody wants to read 200 lines for such a basic question.

Regards

Bruno

 Jessica Soares da Rocha (jessicarocha) said on 2019-05-12: #2

I’m so sorry for the script. I want to know how to calculate the total force acting in a clump. I tried this script but it didn’t work:

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)

 Jérôme Duriez (jduriez) said on 2019-05-13: #3

Hi,

YADE clumps are themselves bodies, hence have an id (let's say it's clumpId), and a resulting force that can be accessed through

O.forces.f(clumpId)

PS : I did not really look at your script either

 Bruno Chareyre (bruno-chareyre) said on 2019-05-13: #4

#2 is not a script.
Please check [1], specifically point 3.