← Back to team overview

yade-users team mailing list archive

[Question #699428]: Cannot simulate periodic box with walls

 

New question #699428 on Yade:
https://answers.launchpad.net/yade/+question/699428

Dear YADE users,

My goal is to simulate a representative volume element under uniaxial compression. I want to compare two cases:
1) a fully periodic simulation in which the particle coordinates are scaled
2) a simulation in which the particles are compressed by walls at the top and bottom of the box, with periodicity in all other directions.

Case 1 is easy to simulate, using tutorial example 06-periodic-triaxial-test.py. To simulate case 2 I used example 06-periodic-triaxial-test.py as a starting point and started integrating elements from 03-oedometric-test.py.

I can insert a wall in the periodic box of example 06. However, as soon as I also add the sphere-wall interaction to the collider, the simulation no longer runs. YADE gives no error messages but refuses to iterate past iteration 1.

The input script I would like to run is (note that although this script is more complicated, I already run into issues when only adding the wall and sphere-wall interactions to example 06):

from __future__ import print_function
sigmaIso = -1e5
sigmaMin = -1e3

#import matplotlib                                                                                                                                                                            
#matplotlib.use('Agg')                                                                                                                                                                        

# generate loose packing                                                                                                                                                                      
from yade import pack, qt, plot

O.periodic = True
sp = pack.SpherePack()
# uniform distribution                                                                                                                                                                        
sp.makeCloud((0, 0, 0), (2, 2, 2), rMean=.1, rRelFuzz=.3, periodic=True)

# setup periodic boundary, insert the packing                                                                                                                                                 
sp.toSimulation()

O.bodies.append(wall(0.0, axis=2, sense=1))

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()], allowBiggerThanPeriod=True),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
                        [Ip2_FrictMat_FrictMat_FrictPhys()],
                        [Law2_ScGeom_FrictPhys_CundallStrack()]),
        NewtonIntegrator(damping=.2),
        PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker')
]
O.dt = .25 * PWaveTimeStep()


def checkUnbalanced():
        if O.iter < 5000:
                return
        if unbalancedForce() > 0.1:
                return
        O.bodies.append(wall(2.0, axis=2,sense=-1))
        global plate
        plate.state.vel = (0,0,-0.01)
        O.engines = O.engines + [PyRunner(command='addPlotData()',iterPeriod=200)]
        checker.command = 'unloadPlate()'

def unloadPlate():
        if abs(O.forces.f(plate.id)[2]/(2*2)) > abs(sigmaIso):
                plate.state.vel *= -1
                checker.command = 'stopUnloading()'

def stopUnloading():
        if abs(O.forces.f(plate.id)[2]/(2*2)) < abs(sigmaMin):
                plot.saveDataTxt(O.tags['d.id']+'.txt')
                O.pause()

def addPlotData():
        if not isinstance(O.bodies[-1].shape,Wall):
                plot.addData()
                return
        Fz = O.forces.f(plate.id)[2]
        plot.addData(
                unbalanced=unbalancedForce(),
                i=O.iter,
                sz=abs(fz)/(2*2),
                ev=abs(plate.state.pos[2]*2*2-2*2*2)/(2*2*2)
        )


# define what to plot                                                                                                                                                                         
plot.plots = {'i': ('unbalanced',),'ev': ('sz'),}
# show the plot                                                                                                                                                                               
plot.plot()

O.run()

With kind regards,
Danny van der Haven



-- 
You received this question notification because your team yade-users is
an answer contact for Yade.