yade-users team mailing list archive

[Antoine FAULCONNIER <question697909@xxxxxxxxxxxxxxxxxxxxx>]

Question #697909 on Yade changed:
https://answers.launchpad.net/yade/+question/697909

    Status: Needs information => Open

Antoine FAULCONNIER gave more information on the question:
Hello Jan Stránský ,

Sorry, I am quite new to yade and this forum. I will add precisions.
The simulation is divided in 3 steps : (i) gravity deposition inside the box (ii) compaction of the grains by translation of the top wall and (iii) shearing of the pack of grains by rotating two of the side walls. 
The execution is a bit long so I save after each step, to change parameters, reload, etc...

Here is my code where the unexpected angle of the wall happends (sorry
it is a bit long) :

# -*- coding: utf-8 -*-

from yade import plot,utils
from yade import pack
import math
import random
import numpy
import os.path

model = {}
model['THETA']     = 0.031  #yade.params.table.theta
model['COMPACITY'] = 0.55  #yade.params.table.compacity

## Define the filename for saved state. By including paramater values in the name
## we make sure that a new state is generated for each parameter set
step1_deposition  = "Step1_deposition.yade.gz"
step2_compression = "Step2_compaction_compacity="+str(model['COMPACITY'])+".yade.gz"
step3_shearing = "Step3_shearing_compacity="+str(model['COMPACITY'])+"theta="+str(model['THETA'])+".yade.gz"
## Check if a saved state exists and proceed to function/variable declarations as usual.
import os.path
savedState1 = os.path.exists(step1_deposition)
savedState2 = os.path.exists(step2_compression)
savedState3 = os.path.exists(step3_shearing)

# general model parameters 
model['f']              = 250                                # [Hz]    excitation frequency     
model['T']              = 1/model['f']                    # [s]     excitation period 
model['omega']          = 2*pi*model['f']   # [rad/s] excitation pulsation 
model['e_n']            = 0.63	                       #         normal viscous damping
model['e_s']            = 0.63                  #1                 #         tangential viscous damping
model['friction_angle'] = 0.197             # friction angle, coulomb coeff \mu = tan(friction_angle)
model['g']              = 9.81
#model['COMPACITY']      = 0.6

# grain's parameters
grains = {}
grains['radius']   = 0.00075       # [m]
grains['young']    = 70.0e9        # [Pa]
grains['poisson']  = 0.22          # [1]
grains['density']  = 2500.0        # [kg/m³]
grains['friction'] = 0.197         # [1]

# wall parameters
box = {}
box['initial_height']   = 0.017        # [m] initial heigth of the box before compaction in z direction
box['length']           = 0.017        # [m] length of the container in x direction
box['width']            = 0.017        # [m] width of the container in y direction

## MATERIAL CREATION   
O.materials.append(FrictMat(young = grains['young'], poisson = grains['poisson'], frictionAngle = grains['friction'], density = grains['density']  , label = 'spmat'))   # sphere material

## sphere packing creation
grains_box_offset = 0.001
pred = pack.inAlignedBox((-box['width']/2+grains_box_offset ,-box['width']/2+grains_box_offset   ,0+grains_box_offset ), (box['width']/2 -grains_box_offset ,box['width']/2 -grains_box_offset  ,box['initial_height']-grains_box_offset ))
O.bodies.append(pack.regularHexa(pred, radius=grains['radius'] , gap=0.0001, material = 'spmat'))

## random distribution of radii
Dev     = 0.00005
media   = 0.00075
val_max = media + Dev
val_min = media - Dev

radii = numpy.linspace(val_min, val_max, num=50, endpoint=True,
retstep=False)

for b in O.bodies:
    if isinstance(b.shape,Sphere):
        b.shape.radius = random.choice(radii)

## total volume and number of the spheres
V_spheres = 0
number = 0

for b in O.bodies:
    if isinstance(b.shape, Sphere):
        V_spheres = V_spheres + 4/3*pi*b.shape.radius**3
        number = number + 1

## height from Compacity
box['final_height'] = V_spheres/model['COMPACITY']*1/(box['length']*box['width'] )  #final height of the box after compression of the grains

## Creation of the container
Down    = O.bodies.append(geom.facetBox((0,0,0),(box['length'],box['width']/2 ,0)))
Right   = O.bodies.append(geom.facetBox((box['length']/2,-box['width']/2 ,box['initial_height']/2), (0,box['width'],box['initial_height']/2)))
Left    = O.bodies.append(geom.facetBox((-box['length']/2,-box['width']/2 ,box['initial_height']/2),(0,box['width'],box['initial_height']/2)))
Front   = O.bodies.append(geom.facetBox((0,-box['width']/2,box['initial_height']),(box['length']/2,0,box['initial_height'])))
Bottom  = O.bodies.append(geom.facetBox((0,box['width']/2,box['initial_height']),(box['length'],0,box['initial_height'])))
Cover   = O.bodies.append(geom.facetBox((0,0,box['initial_height']),(box['length'],box['width']/2 ,0)))

## FUNCTIONS FOR SAVING DATA 
def plotData():
    plot.addData(
    t        = O.time,
    Edf      = Mindlin.frictionDissipation,
    Eds      = Mindlin.shearDampDissip, 
    Edn      = Mindlin.normDampDissip,
    Etot     = O.energy.total(),
    Ed       = Mindlin.frictionDissipation + Mindlin.shearDampDissip + Mindlin.normDampDissip,
    T        = -utils.sumTorques(ids=Right,axis = (0,0,1),axisPt = (box['length']/2,-box['width']/2 ,box['initial_height']/2)),
    TL       = -utils.sumTorques(ids=Left,axis = (0,0,1),axisPt = (-box['length']/2,-box['width']/2,box['initial_height']/2)),
    AngleR    = O.bodies[Right[0]].state.ori[2],
    AngleL    = O.bodies[Left[0]].state.ori[2],
    AngVel   = O.bodies[Right[0]].state.angVel[2],
    Pressure = ((O.forces.f(O.bodies[Cover[0]].id)[2] + O.forces.f(O.bodies[Cover[1]].id)[2])/(0.0002))*1e-6,
    Top_Wall_pos = O.bodies[Cover[0]].state.pos[2],
    **O.energy
    )
    if O.iter > 1000000:
        plot.saveDataTxt('Side'+str(model['THETA'])+'compac'+str(model['COMPACITY'])+'test_save.dat', vars=('t', 'Edf', 'Eds', 'Edn', 'T', 'TL','AngVel','AngleR','AngleR','Pressure','Top_Wall_pos'))

## engine
O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys(frictAngle = model['friction_angle'], en=model['e_n'] , es=1, label='Ip2')],     # damping is present but no friction to reach quickly a balanced state
        [Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=False, calcEnergy=True, label='Mindlin')]
    ),
    NewtonIntegrator(gravity=(0, 0, -model['g']), damping=0.0),
    TranslationEngine(dead = True, translationAxis=[0, 0, 1], velocity=-0.002, ids=Cover, label='pres'),
    HarmonicRotationEngine(dead=True, f =model['f'], A = 2*model['THETA'],ids=Right,rotateAroundZero=True,rotationAxis=(0, 0, 1),zeroPoint=(box['length']/2,-box['width']/2,0),label ='rotR'),    #fi = 1*pi/2,
    HarmonicRotationEngine(dead=True, f =model['f'], A = 2*model['THETA'], ids=Left, rotateAroundZero=True, rotationAxis=(0, 0, 1), zeroPoint=(-box['length']/2, -box['width']/2, 0), label = 'rotL'), #fi = 1*pi/2,
    PyRunner(dead = True, command = 'startPressure()', iterPeriod = 100, label = 'pressure_function'),
    PyRunner(dead = False, command = 'plotData()', iterPeriod = 1000, label = 'data')
]

O.dt = PWaveTimeStep()
flag = True
O.trackEnergy = True

## pressure control fonction
def startPressure():
    pres.dead = False
    pres.velocity = -0.5 
    if (O.bodies[Cover[0]].state.pos[2] - top_sphere)< 0.0001:
        pres.velocity = -0.001
    elif O.bodies[Cover[0]].state.pos[2] <= box['final_height']:
        pres.dead = True       

## simulation process
if not savedState2:            # is there a compressed state ?
    if not savedState1:		# is there an initial deposited state ?
        print("no gravity deposited state found - running script from beginning, starting with gravity deposition")
        O.run(1000000,True)
        while O.iter < 100000 and unbalancedForce() > 0.01:
            O.run(100000,True)
        print("Deposition achieved - save deposition then proceed with compaction")            
        O.save(step1_deposition)
    else:
        print("gravity deposited state found - proceeding with the compaction")
        O.load(step1_deposition)
    ## compression + resting part
    top_sphere = max([b.state.pos[2] for b in O.bodies if isinstance(b.shape, Sphere)]) 
    ## calculatig approximately the number of time steps that are required to place the cover according to the compacity required
    nb_dt_press = round(((box['initial_height']-top_sphere)/0.5+(top_sphere-box['final_height'])/0.001)/O.dt)
    pressure_function.dead = False
    O.run(nb_dt_press+200000,True)
    pressure_function.dead = True
    Ip2.frictAngle = model['friction_angle']
    O.materials[0].frictionAngle = model['friction_angle']
#    O.materials[1].frictionAngle = model['friction_angle']
    Ip2.en = model['e_n']
    Ip2.es = model['e_s']
    print(((O.forces.f(O.bodies[Cover[0]].id)[2] + O.forces.f(O.bodies[Cover[1]].id)[2])/(1))*1)
    pres.dead = True
    O.bodies[Cover[0]].state.blockedDOFs = 'xyzXYZ'
    O.bodies[Cover[0]].state.vel = (0,0,0)
    O.bodies[Cover[1]].state.blockedDOFs = 'xyzXYZ'
    O.bodies[Cover[1]].state.vel = (0,0,0)
    print("compaction and resting phase finished - save compaction then proceed with shearing")        
    O.save(step2_compression)     
else:
    print("compressed state found - will now proceed with shearing") 
    O.load(step2_compression)    
rotR.dead = False
rotL.dead = False
N_cycle = 2            
O.run(N_cycle*round(model['T']/O.dt),True) 
print("Shearing phase achieved - proceed with saving")
O.save(step3_shearing)

Thanks for your help
Antoine

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