yade-users team mailing list archive

[Jessica Soares da Rocha <question679808@xxxxxxxxxxxxxxxxxxxxx>]

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.