# yade-users team mailing list archive

## Re: [Question #697909]: unexpected offset angle with the HarmonicRotationEngine

Question #697909 on Yade changed:

Status: Needs information => Open

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 -*-

import math
import random
import numpy
import os.path

model = {}

## 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
## 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['young']    = 70.0e9        # [Pa]
grains['poisson']  = 0.22          # 
grains['density']  = 2500.0        # [kg/m³]
grains['friction'] = 0.197         # 

# 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 ))

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):

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

for b in O.bodies:
if isinstance(b.shape, Sphere):
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():
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].state.ori,
AngleL    = O.bodies[Left].state.ori,
AngVel   = O.bodies[Right].state.angVel,
Pressure = ((O.forces.f(O.bodies[Cover].id) + O.forces.f(O.bodies[Cover].id))/(0.0002))*1e-6,
Top_Wall_pos = O.bodies[Cover].state.pos,
**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
),
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.velocity = -0.5
if (O.bodies[Cover].state.pos - top_sphere)< 0.0001:
pres.velocity = -0.001
elif O.bodies[Cover].state.pos <= box['final_height']:

## 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")
## compression + resting part
top_sphere = max([b.state.pos 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)
O.run(nb_dt_press+200000,True)
Ip2.frictAngle = model['friction_angle']
O.materials.frictionAngle = model['friction_angle']
#    O.materials.frictionAngle = model['friction_angle']
Ip2.en = model['e_n']
Ip2.es = model['e_s']
print(((O.forces.f(O.bodies[Cover].id) + O.forces.f(O.bodies[Cover].id))/(1))*1)
O.bodies[Cover].state.blockedDOFs = 'xyzXYZ'
O.bodies[Cover].state.vel = (0,0,0)
O.bodies[Cover].state.blockedDOFs = 'xyzXYZ'
O.bodies[Cover].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")
N_cycle = 2
O.run(N_cycle*round(model['T']/O.dt),True)
print("Shearing phase achieved - proceed with saving")
O.save(step3_shearing)