ConnectionFinder problem and Bonding broke after the first frame

Asked by Wasin Meesuay on 2020-11-01

According to my previous question (https://answers.launchpad.net/esys-particle/+question/693648), I found out that using the ConnectionFinder to define bonding between particles need to set the MaxDist up to the Minimum radius to make all particles bonded to each other.

> I am very curious that, normally the MaxDist in ConnectionFinder usually had low value (around 0.005 - 1) and can bond all particles together, but in my case it seem to need a large value (around Rmin which is up to 10-40). Why the MaxDist is this large?

Then, after I define all the bond and conducted a 2D biaxial testing, all of bonding just appear properly in the first frame, but in the next frame most of them just disappear (Only a few just left and also broke in next 2-3 step).

First, I think it might be the same with this problem (https://answers.launchpad.net/esys-particle/+question/245754) , so I tried to

> (1.) add linear and rotational damping
> (2.) reduce the timestep size 100 lower than the critical size
> (3.) remove the confining pressure in the beginning

but the bond still broke after the first frame. The only way that I can do to make the bond broke the least is to reduce the Rmin of the particle. However, the number of broken bond still more than half in the second frame without applying any confining pressure.

How do I fix this problem? Do I triggered any bugs?

Thank you for an advance .
Wasin

Here is the script of my biaxial test simulation
####################################

#import the appropriate ESyS-Particle modules:
from esys.lsm import *
from esys.lsm.util import *
from esys.lsm.geometry import *
from WallLoader import WallLoaderRunnable
from ServoWallLoader import ServoWallLoaderRunnable

#Set Parameters
Rmin = 40.0
Rmax = 200.0
Xlength = 4000
Ylength = 3000
YoungMod = 2.00e+9
Poisson = 0.30
Cohes = 20.0e+6
TanPhi = 1.00
MuSta = 0.60
MuDyn = 0.40
Rho = 2500.00
ConfiningP = 10.0e+6

TotalTimestep = 2500000

#instantiate a simulation object:
sim = LsmMpi (numWorkerProcesses = 3, mpiDimList = [3,1,1])

#initialise the neighbour search algorithm:
sim.initNeighbourSearch (
   particleType = "RotSphere",
   gridSpacing = 2.2 * Rmax,
   verletDist = Rmin/5.0
)

#set the number of timesteps and timestep increment:
sim.setNumTimeSteps (TotalTimestep)
sim.setTimeStepSize (1.0e-04)

#enforce two-dimensional computations:
sim.force2dComputations (True)

#specify the spatial domain for the simulation
domain = BoundingBox(Vec3(-0.5 * Xlength,-0.5 * Ylength,0), Vec3( 1.5 * Xlength , 1.5 * Ylength, 0))
sim.setSpatialDomain (domain)

#create a prism of spherical particles:
geoRandomBlock = RandomBoxPacker (
   minRadius = Rmin,
   maxRadius = Rmax,
   cubicPackRadius = 2.2000,
   maxInsertFails = 50000,
   bBox = BoundingBox(
      Vec3(0.0, 0.0, 0.0),
      Vec3(Xlength, Ylength, 0.0)
   ),
   circDimList = [False, False, False],
   tolerance = 1.0000e-05
)
geoRandomBlock.generate()
geoRandomBlock_particles = geoRandomBlock.getSimpleSphereCollection()

#add the particles to the simulation object:
for pp in geoRandomBlock_particles:
    centre = pp.getPosn()
    radius = pp.getRadius()
    Y = centre[1]
    if (Y < Rmin/2): # particle is near the base (tag=2)
        pp.setTag (2)
    elif (Y > (Ylength-Rmax/2)): # particle is near the top (tag=3)
        pp.setTag (3)
    else: # particle is inside the shear cell (tag=1)
        pp.setTag (1)
    sim.createParticle(pp)

#set the density of all particles:
sim.setParticleDensity (
   tag = 1,
   mask = -1,
   Density = Rho
)

#set the density of all particles:
sim.setParticleDensity (
   tag = 2,
   mask = -1,
   Density = Rho
)

#set the density of all particles:
sim.setParticleDensity (
   tag = 3,
   mask = -1,
   Density = Rho
)

#bond particles together with bondTag = 1:
sim.createConnections(
   ConnectionFinder(
      maxDist = Rmin,
      bondTag = 1,
      pList = geoRandomBlock_particles
   )
)

#create a wall at the bottom of the model:
sim.createWall (
   name = "bottom_wall",
   posn = Vec3(0.0, 0.0, 0.0),
   normal = Vec3(0.0, 1.0, 0.0)
)

#create a wall at the top of the model:
sim.createWall (
   name = "top_wall",
   posn = Vec3(0.0, Ylength, 0.0),
   normal = Vec3(0.00, -1.0, 0.0)
)

#create a wall at the bottom of the model:
sim.createWall (
   name = "left_wall",
   posn = Vec3(0.0, 0.0, 0.0),
   normal = Vec3(1.0, 0.0, 0.0)
)

#create a wall at the bottom of the model:
sim.createWall (
   name = "right_wall",
   posn = Vec3(Xlength, 0.0, 0.0),
   normal = Vec3(-1.0, 0.0, 0.0)
)

#create rotational elastic-brittle bonds between particles:
pp_bonds = sim.createInteractionGroup (
   BrittleBeamPrms(
      name="pp_bonds",
      youngsModulus= YoungMod,
      poissonsRatio= Poisson,
      cohesion= Cohes,
      tanAngle= TanPhi,
      tag=1
   )
)

#initialise frictional interactions for unbonded particles:
sim.createInteractionGroup (
   FrictionPrms(
      name="friction",
      youngsModulus= YoungMod,
      poissonsRatio= Poisson,
      dynamicMu= MuDyn,
      staticMu= MuSta
   )
)

#create an exclusion between bonded and frictional interactions:
sim.createExclusion (
   interactionName1 = "pp_bonds",
   interactionName2 = "friction"
)

#specify elastic repulsion from the bottom wall:
sim.createInteractionGroup (
   NRotBondedWallPrms (
      name = "bottom_wall_repel",
      wallName = "bottom_wall",
      normalK = YoungMod,
      particleTag = 2
   )
)

#specify elastic repulsion from the top wall:
sim.createInteractionGroup (
   NRotBondedWallPrms (
      name = "top_wall_repel",
      wallName = "top_wall",
      normalK = YoungMod,
      particleTag = 3
   )
)

#specify elastic repulsion from the top wall:
sim.createInteractionGroup (
   NRotElasticWallPrms (
      name = "left_wall_repel",
      wallName = "left_wall",
      normalK = YoungMod
   )
)

#specify elastic repulsion from the top wall:
sim.createInteractionGroup (
   NRotElasticWallPrms (
      name = "right_wall_repel",
      wallName = "right_wall",
      normalK = YoungMod
   )
)

#add translational viscous damping:
sim.createInteractionGroup (
   LinDampingPrms(
      name="damping1",
      viscosity=2,
      maxIterations=100
   )
)

#add rotational viscous damping:
sim.createInteractionGroup (
   RotDampingPrms(
      name="damping2",
      viscosity=2,
      maxIterations=100
   )
)

#add a wall loader to move the top wall:
wall_loader1 = WallLoaderRunnable(
   LsmMpi = sim,
   wallName = "right_wall",
   vPlate = Vec3 (-0.1, 0.0, 0.0),
   startTime = TotalTimestep/5,
   rampTime = TotalTimestep/4
)
sim.addPreTimeStepRunnable (wall_loader1)

#add a wall loader to move the bottom wall:
wall_loader2 = WallLoaderRunnable(
   LsmMpi = sim,
   wallName = "left_wall",
   vPlate = Vec3 (0.1, 0.0, 0.0),
   startTime = TotalTimestep/5,
   rampTime = TotalTimestep/4
)
sim.addPreTimeStepRunnable (wall_loader2)

#add ServoWallLoaderRunnables to apply constant normal stress:
servo_loader1 = ServoWallLoaderRunnable(
   LsmMpi = sim,
   interactionName = "top_wall_repel",
   force = Vec3 (0.0,(-0.5*Xlength*ConfiningP), 0.0),
   startTime = 0,
   rampTime = 400000
)
sim.addPreTimeStepRunnable (servo_loader1)

#add ServoWallLoaderRunnables to apply constant normal stress:
servo_loader2 = ServoWallLoaderRunnable(
   LsmMpi = sim,
   interactionName = "bottom_wall_repel",
   force = Vec3 (0.0, (0.5*Xlength*ConfiningP), 0.0),
   startTime = 0,
   rampTime = 400000
)
sim.addPreTimeStepRunnable (servo_loader2)

#create a FieldSaver to store number of bonds:
sim.createFieldSaver (
   InteractionScalarFieldSaverPrms(
      interactionName="pp_bonds",
      fieldName="count",
      fileName="nbonds.dat",
      fileFormat="SUM",
      beginTimeStep=0,
      endTimeStep=TotalTimestep,
      timeStepIncr=1000
   )
)

#create a FieldSaver to store the total kinetic energy of the particles:
sim.createFieldSaver (
   ParticleScalarFieldSaverPrms(
      fieldName="e_kin",
      fileName="ekin.dat",
      fileFormat="SUM",
      beginTimeStep=0,
      endTimeStep=TotalTimestep,
      timeStepIncr=1000
   )
)

#create a FieldSaver to store potential energy stored in bonds:
sim.createFieldSaver (
   InteractionScalarFieldSaverPrms(
      interactionName="pp_bonds",
      fieldName="potential_energy",
      fileName="epot.dat",
      fileFormat="SUM",
      beginTimeStep=0,
      endTimeStep=TotalTimestep,
      timeStepIncr=1000
   )
)

#create a FieldSaver to wall positions:
posn_saver = WallVectorFieldSaverPrms(
   wallName=["bottom_wall", "top_wall", "left_wall", "right_wall"],
   fieldName="Position",
   fileName="out_Position.dat",
   fileFormat="RAW_SERIES",
   beginTimeStep=0,
   endTimeStep=TotalTimestep,
   timeStepIncr=1000
)
sim.createFieldSaver(posn_saver)

#create a FieldSaver to wall forces:
force_saver = WallVectorFieldSaverPrms(
   wallName=["bottom_wall", "top_wall", "left_wall", "right_wall"],
   fieldName="Force",
   fileName="out_Force.dat",
   fileFormat="RAW_SERIES",
   beginTimeStep=0,
   endTimeStep= TotalTimestep,
   timeStepIncr=1000
)
sim.createFieldSaver(force_saver)

#add a CheckPointer to store simulation data:
sim.createCheckPointer (
   CheckPointPrms (
      fileNamePrefix = "snapshot",
      beginTimeStep = 0,
      endTimeStep = TotalTimestep,
      timeStepIncr = 100000
   )
)

#execute the simulation:
sim.run()

Question information

Language:
English Edit question
Status:
Solved
For:
ESyS-Particle Edit question
Assignee:
No assignee Edit question
Solved by:
Dion Weatherley
Solved:
2020-12-03
Last query:
2020-12-03
Last reply:
2020-11-09
Wasin Meesuay (kakkunno) said : #1

I also pop up an idea that might solve the broken bond problem by creating a new bonding after the simulation start for a while.
Is there any chance that I could create the bond during the simulation??
If possible, it would be grateful that you could guide me through some pieces of code. (I quite very newbie on esys-particle)

Thank you for an advance krub :) ( ^__/\__ ^)

Best Dion Weatherley (d-weatherley) said : #2

Hi Wasin,

I think the problem is arising due to RandomBoxPacker and ConnectionFinder. Try changing your definition of RandomBoxPacker to the following:

#create a prism of spherical particles:
geoRandomBlock = RandomBoxPacker (
   minRadius = Rmin,
   maxRadius = Rmax,
   cubicPackRadius = 2.2000*Rmax,
   maxInsertFails = 50000,
   bBox = BoundingBox(
      Vec3(0.0, 0.0, 0.0),
      Vec3(Xlength, Ylength, 0.0)
   ),
   circDimList = [False, False, False],
   tolerance = 1.0000e-05
)
geoRandomBlock.generate()
geoRandomBlock_particles = geoRandomBlock.getSimpleSphereCollection()

Notice that I have changed cubicPackRadius to =2.2*Rmax, as opposed to 2.2.

Also, change the ConnectionFinder to:
#bond particles together with bondTag = 1:
sim.createConnections(
   ConnectionFinder(
      maxDist = 0.1,
      bondTag = 1,
      pList = geoRandomBlock_particles
   )
)

Here I have reduced the maxDist to 0.1 from Rmin.

When I run this example, I get around 450 particles and 800 bonds. One bond breaks in the first timestep but otherwise the simulation appears to be stable. Perhaps reduce maxDist to 0.05 or smaller to avoid this single bond breakage?

Cheers,

Dion

Wasin Meesuay (kakkunno) said : #3

Thanks Dion Weatherley, that solved my question.