Yade

Stress remains 0 during the simulation of a direct shear test

Asked by Vinicius de Souza Lima

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)

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Karol Brzezinski (kbrzezinski) said : 		#1

Hi,

I noticed you already asked this question in Gitlab (which is now recommended). So, please note that I posted the full answer only there.

https://gitlab.com/yade-dev/answers/-/issues/3876

Cheers,
Karol

Can you help with this problem?

Provide an answer of your own, or ask Vinicius de Souza Lima for more information if necessary.

To post a message you must log in.