yade-users team mailing list archive

Thread
Date
[Question #680746]: How to plot the total Force on a Clump

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: Jessica Soares da Rocha <question680746@xxxxxxxxxxxxxxxxxxxxx>
Date: Fri, 10 May 2019 22:17:29 -0000
Reply-to: question680746@xxxxxxxxxxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx
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.