unexpected offset angle with the HarmonicRotationEngine

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

[1] https://www.yade-dem.org/wiki/Howtoask

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

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

> I am quite new to yade and this forum

welcome :-)
and please read [1] carefully, e.g. also add what is your Yade version (sometimes there are changes, that might be important for the problem) etc.

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

[1] https://www.yade-dem.org/wiki/Howtoask

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

Cheers
Antoine

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