apply rigid body rotation of a periodic cell

Hi Ning,

>
> I want to apply a rigid body rotation on a periodic cell meantime maintain
> the configuration (relative particle position and contact force unchanged
> but only rotated). If I use O.cell.velGrad=Matrix3(0,sin(theta),0,
> -sin(theta),0,0, 0,0,0), the cell rotates but the configuration are also
> changed. Is there a way currently available to do the pure rotation?
>
>
the problem here is that the matrix you udefined does not produce pure
rotation. The change of configuration of cell (trsf) is expressed [1]:

trsfNew = trsfOld + velGrad*dt*trsfOld

so if you know how trsf should chenge during the time (rotate in your
case), you can define (after some matrix algebra) propper velGrad.

cheers
Jan

[1]
http://bazaar.launchpad.net/~yade-pkg/yade/git-trunk/view/head:/core/Cell.cpp#L9

I think velGrad=Matrix3(0,sin(theta),0, -sin(theta),0,0, 0,0,0) is correct asymptoticaly for small spin*dt.
If you really need a large rotation in one step, I'm afraid it is not possible currently (Jan's trick sounds good, except that the particles will not be rotated - it could be implemented).

Hi guys, for the rotation of particles with PBC, I've done the homework for you: http://woodem.eu/doc/theory/leapfrog.html#motion-in-uniformly-deforming-space . Cheers, vaclav

Vaclav,
Good to know, but... I think your algorithm will not work combined with Jan's idea #1 because it still assumes a velocity gradient incremented over many small timesteps (this is also a requirement for equating rotation and cross product).
Your spin definition is only ok with #2 (are you sure about the W/2, why not just W?).

@Bruno: OK, I get what you mean.

EH, not sure about W/2 now, I guess I just got it from Truesdell&Noll or such, but don't remember to be honest. I can only see that it works as it should (spheres are rotated with the space) but there might be other 1/2 factor somewhere which is missing or such... :| tricky business.

If someone knows how to smuggle spin in to the aspherical integrator, I will be glad to know.

Hi Ning,

could you describe your physical situation in which you want to apply a rigid body rotation on a periodic cell?

The idea is quite interesting, however in general case this is not possible, because fictitious forces depend on distance from the axis of rotation. So in every periodic cell the forces will be different (they have different distance from the axis of rotation, and you can have this axis in one cell only). So suddenly due to rotation each cell in an infinite arrangement of cells is different.

Check http://en.wikipedia.org/wiki/Fictitious_force#Rotating_coordinate_systems

best regards
Janek Kozicki

Thank you for the kind replies,

I'm using YADE for coupled FEM/DEM simulation where DEM packings serve as RVEs in Gauss points of FEM mesh. Some papers describe this idea [1,2,3]. Now I'm dealing with a problem with cylindrical geometry. I want the initial condition to be axial-symmetric. So I'd to rotate the initial RVE packings according to their positions in the problem domain because the packings are not perfectly isotropic (only ~400 grains in each packing).

But later I think I can solve the problem (FEM level) in the cylindrical coordinate. Then I may not need to rotate the RVE packings. So currently I'm working on the FEM formulation in the cylindrical coordinate. And sorry for the late reply.

[1] http://link.springer.com/article/10.1007%2Fs10035-011-0255-6
[2] http://dx.doi.org/10.1063/1.4812158
[3] http://dx.doi.org/10.1063/1.4812151

Best regards,
Ning

Hi Ning,
what I understood is that what you actually want is to place you packing to
already rotated cell. This is possible with SpherePack.toSimulation
function [1]
cheers
Jan

[1]
https://yade-dem.org/doc/yade.pack.html#yade.pack.SpherePack_toSimulation

2014-05-06 6:06 GMT+02:00 ceguo <email address hidden>:

> Question #247806 on Yade changed:
> https://answers.launchpad.net/yade/+question/247806
>
> ceguo posted a new comment:
> Thank you for the kind replies,
>
> I'm using YADE for coupled FEM/DEM simulation where DEM packings serve
> as RVEs in Gauss points of FEM mesh. Some papers describe this idea
> [1,2,3]. Now I'm dealing with a problem with cylindrical geometry. I
> want the initial condition to be axial-symmetric. So I'd to rotate the
> initial RVE packings according to their positions in the problem domain
> because the packings are not perfectly isotropic (only ~400 grains in
> each packing).
>
> But later I think I can solve the problem (FEM level) in the cylindrical
> coordinate. Then I may not need to rotate the RVE packings. So currently
> I'm working on the FEM formulation in the cylindrical coordinate. And
> sorry for the late reply.
>
> [1] http://link.springer.com/article/10.1007%2Fs10035-011-0255-6
> [2] http://dx.doi.org/10.1063/1.4812158
> [3] http://dx.doi.org/10.1063/1.4812151
>
> Best regards,
> Ning
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi Jan,

Actually I need all the RVE packings possess exactly same initial condition except they are rigidly rotated. Because I'm trying to prepare an axi-symmetric problem domain. See the illustration figure:

https://dl.dropboxusercontent.com/u/7120187/illust.png

cheers,
Ning

@Ning,
I see two possible simple solutions to your problem:
One needs cpu time (see example below):
dtheta=0.01 #something small (aka rotation per timestep)
thetaDot=dtheta/O.dt
velGrad=Matrix3(0,thetaDot,0, -thetaDot,0,0, 0,0,0)
N=theta/dtheta
O.run(N) #now everything is rotated by theta

The other is instantaneaous: instead of rotating the packing, rotate the imposed velocity gradient in order to define it in the cells reference frame. For theta=45°:
q=Quaternion().setFromTwoVectors((1,0,0),(1,1,0))
velGrad=Matrix3(1,0,0,0,-1,0,0,0,0) #pure shear, written in the fixed cartesian frame
q.toRotationMatrix()*velGrad*q.inverse.toRotationMatrix()
Matrix3(-2.220446049250313e-16,1,0, 1,2.220446049250313e-16,0, 0,0,0) #pure shear, in the cells reference frame

Actually, you may not have to rotate anything if the strain from FEM is defined in an axisymmetric frame, since r,theta are always aligned on the periodic cell axis...

@Vaclav
My impression is that W/2 comes out in simple shear (and so it is visible in many books - and correct), as the skew-symmetric part of L when L is:
0 W
0 0
However, when the gradient is antisymmetric in the first place: 1/2(L+LT) = L, with L=
 0 W
-W 0
The 1/2 disappears, and the spin is just W.

Whatever it is, it will work when you look at a simulation. Even with yade, without defining any spin, they rotate correctly in order to maintain equilibrium (exception for a few spheres without contacts, which clearly do not rotate correctly).

Example of rotating (modified from test/periodic-simple-shear.py):
============

from yade import pack,qt
O.periodic=True
O.cell.hSize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
radius=5e-3
sp=pack.SpherePack().makeCloud((0,0,0),(.1,.1,.1),5e-3,.2,500,periodic=True)

num=sp.makeCloud((0,0,0),(.1,.1,.1),radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
O.bodies.append([sphere(s[0],s[1]) for s in sp])
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=.05*radius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ), PeriTriaxController(dynCell=True,mass=0.2,maxUnbalanced=0.01,relStressTol=0.02,goal=[-1e4,-1e4,0],stressMask=3,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='triaxDone()',label='triax'),
 NewtonIntegrator(damping=.2),
]

phase=0
def triaxDone():
 global phase
 if phase==0:
  print 'Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff
  print 'Now rotating.'
  O.cell.velGrad=Matrix3(0,1,0, -1,0,0, 0,0,0)
  triax.dead=1
  phase+=1

O.dt=PWaveTimeStep()
Gl1_Sphere.stripes=1

@Bruno: thanks for the spin info. The rotation I meant is about rotating free-standing spheres with no forces applied, like this: http://youtu.be/bzDo4aHIYPE . When the packing is dense and there is friction, then position update will make particles rotate even without any special consideration for spin, if that's what you meant. Cheers, v.

@Vaclav. The movie prooves that the implementation is correct, indeed.
It must be just something I don't understand in the documentation you linked.