# unexpected offset angle with the HarmonicRotationEngine

Hello everyone,

I am shearing a pack of grains that are contained in a parallelepipedic box by imposing a periodic rotation on two of the side walls, the two rotating walls are facing each other, pure shearing.
For that I am using the partial engine "HarmonicRotationEngine" to rotate the walls.
So the rotation angle is sinusoidal in time.

The problem I encounter is that this oscillation of the wall is not centered on zero, there is an offset value of the angle.
I tried to change the value of the phase "fi" in the partial engine, from 0 to pi/2, it was better, still not centered.
I also read in the documentation that by default (no "fi" inputted) the body should rotate around its initial position, which is not the case in my simulation.

Why doesn't it rotate around the initiale position ?

I hope it is clear enough.

Thanks

Antoine

PS : here is the command line for the partial engine:
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')

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Antoine FAULCONNIER
Solved:
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2021-07-09: #1

Hello,

> PS : here is the command line for the partial engine:

please provide a complete MWE [1]
There might be maany reasons, without the actual code we can merely guess.

Cheers
Jan

 Revision history for this message Antoine FAULCONNIER (afaulconnier) said on 2021-07-09 #2

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...
Also I record the angle of the wall in the function plotData() with the commande : AngleR = O.bodies[Right[0]].state.ori[2],

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

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[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
),
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[0]].state.pos[2] - top_sphere)< 0.0001:
pres.velocity = -0.001
elif O.bodies[Cover[0]].state.pos[2] <= 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[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)
O.run(nb_dt_press+200000,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)
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")
N_cycle = 2
O.run(N_cycle*round(model['T']/O.dt),True)
print("Shearing phase achieved - proceed with saving")
O.save(step3_shearing)

Antoine

 Revision history for this message Jan Stránský (honzik) said on 2021-07-09: #3

> I am quite new to yade and this forum

welcome :-)

thanks for the code.

> The execution is a bit long ...
> sorry it is a bit long ...

therefore the code posted should be a MWE [1], M = minimal.
Something like this, a short code focusing just on the problem.
###
bot = (
sphere((0,0,0),1),
sphere((2,0,0),1),
sphere((0,2,0),1),
sphere((2,2,0),1),
)
top = (
sphere((0,0,3),1),
sphere((2,0,3),.5),
sphere((0,2,3),.5),
sphere((2,2,3),1),
)
botIds = O.bodies.append(bot)
topIds = O.bodies.append(top)

harmRot = HarmonicRotationEngine(
ids=topIds,
rotateAroundZero=True,
rotationAxis=(0,0,1),
zeroPoint=(1,1,0),
A=1,
f=1,
)

fResetter,collider,iLoop,timeStepper,newton = O.engines
O.engines = [fResetter,collider,iLoop,harmRot,newton]
O.dt = 1e-6
###

if this is what you are expecting and your code does not work as expected, you can think what is different, add relevant parts to the MWE and iterate similarly.
This way, there is also higher probability that you find the problem yourself.
So please try to create a MWE, i.e. deleting the deposition and compaction and just add a few "representative" hardcoded spheres (i.e. if particles are far away, make them far away, if they are touching, make them touching etc.), deleting not important PyRunners etc.

Cheers
Jan

 Revision history for this message Antoine FAULCONNIER (afaulconnier) said on 2021-07-14: #4

Thanks Jan for your replies, I will try to isolate the problem and come back here after that.

Cheers
Antoine

 Revision history for this message Antoine FAULCONNIER (afaulconnier) said on 2021-07-20 #5

Hello Jan and everyone,

Thank you for the welcoming.

I did not have a lot of time these days but I finally solved my problem, it seems that, the offset amplitude of the sine function guiding the HarmonicRotationEngine is dependant on the phase.

And the origin of the phase is not the moment I start the HarmonicRotationEngine with the command rotengine.dead = False, but the starting moment of the whole simulation. My bad, I did not know.

Anyway I just shifted the phase of the rotation engine so the sine has a phase that is a multiple of pi at the moment I start the rotation engine. And it works fine.

I hope it was understandable.

Thanks for the Help
Antoine