yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #30012
[Question #707913]: Stress remains 0 during the simulation of a direct shear test
New question #707913 on Yade:
https://answers.launchpad.net/yade/+question/707913
Hi! In this simulation, the stress of the top plate(u1) is not increasing in the consolidation phase(before the shearing) although a velocity is imposed. Inside the shear box there is a compact parking of spheres. Can you please tell me if I am doing something wrong?
# -*- coding: utf-8 -*-
from yade import ymport, utils, plot,pack
################## Parâmetros
## copy spheres from the packing into the scene
D_CILINDRO = 59.92e-3
H_CILINDRO = 26.8e-3
X_CUBO = D_CILINDRO
Y_CUBO = H_CILINDRO
Z_CUBO = D_CILINDRO
# Parâmetros das partículas
RAIO = 8e-4
VAR_RAIO = 10e-2 # 10% do Raio
ESFERAS_CELULA = 2000
COR = (0,0,1)
CORPO_DE_PROVA = 'cubo'
PACKING = 'cilindro_0.0009_0.1_2000.sqlite'
DAMPING = 0.5
dtCOEFF = 0.5
normalSTRESS = 50e3 # Pa
normalVEL = 1e-3
shearVEL = 10e-5*X_CUBO # try different values to ensure quasi-static conditions
intR=1#1.263 for X1Y1Z1_5k
DENS=1600
YOUNG=40e6
FRICT=40
ALPHA=0.3
TENS=45e5
COH=45e6
iterMAX=40000
Output='shearTest_'+PACKING+'_K10_D3000_Y50e9A03F18T45e5C45e6_500kPa_shearVel1'
###################### Import of the sphere assembly
### material definition
def sphereMat():
return CohFrictMat(
young = YOUNG,
frictionAngle = radians(FRICT),
normalCohesion = 1,
shearCohesion = 1,
density = DENS,
fragile = False)
def wallMat():
return CohFrictMat(
young = YOUNG,
frictionAngle = 0,
density = DENS,
fragile = False)
pred = pack.inAlignedBox((0,0,0),(X_CUBO,Y_CUBO, Z_CUBO))
#O.cell.hSize = Matrix3(D_CILINDRO, 0, 0,
# 0,D_CILINDRO, 0,
# 0, 0, H_CILINDRO)
sp = SpherePack()
memoizeDb = f'{CORPO_DE_PROVA}_{RAIO}_{VAR_RAIO}_{ESFERAS_CELULA}.sqlite'
sp = pack.randomDensePack(pred,
spheresInCell=2000,
radius=RAIO,
rRelFuzz = VAR_RAIO,
memoizeDb = memoizeDb,
returnSpherePack=True)
spheres = sp.toSimulation(color=(0, 0, 1), material = sphereMat)
## preprocessing to get dimensions
dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf
## initial surface
S0=X*Z
## spheres factory
R=0
Rmax=0
numSpheres=0.
for o in O.bodies:
if isinstance(o.shape,Sphere):
o.shape.color=(0.7,0.5,0.3)
numSpheres+=1
R+=o.shape.radius
if o.shape.radius>Rmax:
Rmax=o.shape.radius
Rmean=R/numSpheres
###################### Criação da caixa cisalhante
thickness=Y/100
oversizeFactor=1.3
### Placas de carregamento
### loading platens
O.bodies.append(utils.box(center=(xinf+X/2.,yinf-thickness+Rmean/10.,zinf+Z/2.),
extents=(oversizeFactor*X/2,thickness,oversizeFactor*Z/2),
material=wallMat,fixed=True))
l1 =O.bodies[-1]
O.bodies.append(utils.box(center=(xinf+X/2.,ysup+thickness-Rmean/10.,zinf+Z/2.),
extents=(oversizeFactor*X/2,thickness,oversizeFactor*Z/2),
material=wallMat,fixed=True))
u1=O.bodies[-1]
### add top box
O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,3*(ysup/4)+2*Rmean,zsup/2),
extents=(thickness,Y/4,oversizeFactor*Z/2),
material=wallMat,
fixed=True))
u2=O.bodies[-1]
O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,3*(ysup/4)+2*Rmean,zsup/2),
extents=(thickness,Y/4,oversizeFactor*Z/2),
material=wallMat,
fixed=True))
u3=O.bodies[-1]
O.bodies.append(utils.box(center=(xsup/2,3*(ysup/4)+2*Rmean,zinf-thickness+Rmean/10),
extents=(oversizeFactor*X/2,Y/4,thickness),
material=wallMat,
fixed=True,
wire=True))
u4=O.bodies[-1]
O.bodies.append(utils.box(center=(xsup/2,3*(ysup/4)+2*Rmean,zsup+thickness-Rmean/10),
extents=(oversizeFactor*X/2,Y/4,thickness),
material=wallMat,
fixed=True,
wire=True))
u5=O.bodies[-1]
### add bottom box
O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,(ysup/4)-2*Rmean,zsup/2),
extents=(thickness,Y/4,oversizeFactor*Z/2),
material=wallMat,
fixed=True))
l2=O.bodies[-1]
O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,(ysup/4)-2*Rmean,zsup/2),
extents=(thickness,Y/4,oversizeFactor*Z/2),
material=wallMat,
fixed=True))
l3=O.bodies[-1]
O.bodies.append(utils.box(center=(xsup/2,(ysup/4)-2*Rmean,zinf-thickness+Rmean/10),
extents=(oversizeFactor*X/2,Y/4,thickness),
material=wallMat,
fixed=True,
wire=True))
l4=O.bodies[-1]
O.bodies.append(utils.box(center=(xsup/2,(ysup/4)-2*Rmean,zsup+thickness-Rmean/10),
extents=(oversizeFactor*X/2,Y/4,thickness),
material=wallMat,
fixed=True,
wire=True))
l5=O.bodies[-1]
###################### Engines
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='aabb'),
Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=intR,
label='ss2d3dg'),
Ig2_Box_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='interactionPhys')],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(label='interactionLaw')]
),
GlobalStiffnessTimeStepper(defaultDt=0.1*PWaveTimeStep(),
timestepSafetyCoefficient=dtCOEFF),
TranslationEngine(ids=[u1.id],
translationAxis=(0.,-1.,0.),
velocity=0.,
label='u1_yTranslation'),
TranslationEngine(ids=[u1.id,u2.id,u3.id,u4.id,u5.id],
translationAxis=(-1.,0.,0.),
velocity=0.,
label='u_xTranslation'),
TranslationEngine(ids=[l1.id,l2.id,l3.id,l4.id,l5.id],
translationAxis=(1.,0.,0.),
velocity=0.,
label='l_xTranslation'),
PyRunner(iterPeriod=1,
initRun=True,
command='servoController()',
label='servo'),
NewtonIntegrator(damping=DAMPING,
gravity=(0,-9.81,0),
label='damper'),
PyRunner(iterPeriod=int(iterMAX/1000),
initRun=True,
command='dataCollector()'),
]
###################### Engines definition ( servoController() and dataCollector() )
shearing=False
sigmaN=0
tau=0
Fs1=0
Fs2=0
Xdispl=0
px0=0
Ydispl=0
py0=l1.state.pos[1]
px0=l2.state.pos[0]
py0=l1.state.pos[1]
prevTranslation=0
n=0
F0 = normalSTRESS
def servoController():
global px0, py0, sigmaN, n, Fn1, Fn2, shearing,u1, u2,u3,l1, l2, l3,Xdispl, Ydispl, tau,F0
sigmaN= abs(O.forces.f(u1.id)[1])/S0 # necessary?
Xdispl = abs(l2.state.pos[0]-px0)
tau = 0.0
placas = [u1,u2,u3,u4,u5,l1,l2,l3,l4,l5]
for placa in placas:
tau += abs(O.forces.f(placa.id)[0])
tau = tau/(2*S0)
if abs(sigmaN-normalSTRESS)/normalSTRESS > 0.01:
if sigmaN < normalSTRESS:
u1_yTranslation.velocity += normalVEL/10
if sigmaN>(normalSTRESS):
u1_yTranslation.velocity -= normalVEL/10
if abs((normalSTRESS-sigmaN)/normalSTRESS) < 0.01 and unbalancedForce() < 15/100:
shearing = True
print(f"sigmaN (Pa) = {sigmaN:.2f};",
f"dif percentual (%) = {(normalSTRESS-sigmaN)/normalSTRESS*100.0:.2f};"
f"velocidade da placa (m/s) = {u1.state.vel[1]:.3f};",
f"unbalancedForce (%) = {unbalancedForce()*100.0:.1f};",
f"tau (Pa) = {tau:.5f}")
if shearing:
u_xTranslation.velocity = shearVEL
l_xTranslation.velocity = shearVEL
def dataCollector():
global Xdispl, Ydispl, tau,sigmaN
Xdispl = abs(l2.state.pos[0]-px0)
Ydispl = abs(l1.state.pos[1]-py0)
yade.plot.addData({'t':O.time,
'i':O.iter,
'Xdispl':Xdispl,
'sigmaN':sigmaN,
'tau':tau,
'Ydispl':Ydispl,
'unbF':utils.unbalancedForce()})
plot.saveDataTxt(Output)
# defines window plot
plot.plots={'i':('sigmaN','tau')}
plot.plot()
################# graphical intervace
from yade import qt
qt.Controller()
qt.View()
v=qt.Renderer()
v.dispScale=(1,1,1) # displacements are scaled for visualization purpose
################# to manage interaction detection factor during the first timestep
O.step()
################# initializes the interaction detection factor to its default value (new contacts will be created between strictly contacting particles only)
ss2d3dg.interactionDetectionFactor=-1.
aabb.aabbEnlargeFactor=-1.
################# simulation starts here
O.run(iterMAX)
--
You received this question notification because your team yade-users is
an answer contact for Yade.
