Yade

How to move O.bodies.append(facet()) to a circle or complex trajectory

Asked by 内山康太郎

Thanks to Jean, I was able to come this far. First of all, thank you Jean.

I am doing an analysis on a direct(simple) shear test.
I have created a program like below.

Previously (https://answers.launchpad.net/yade/+question/706234), I used TranslationEngine(translationAxis()) to make O.bodies.append(facet()) a rectangular trajectory.

This time, I want to use TranslationEngine(translationAxis()) to make O.bodies.append(facet()) become a "circular" or "complex trajectory".

So far, for rectangles, TranslationEngine(translationAxis())'s translationAxis is done at a constant value and only stops when it reaches the specified displacement.
However, with a circular shape, the "translationAxis()" is no longer constant because the "translationAxis()" is made to be a circular trajectory. Also, the (X,Y) coordinates change for each step of the circular shape, so we have to think very hard about it.

So, is there any way to move O.bodies.append(facet()) circularly?
By #circular locus, I mean a locus with respect to a (constant) radius.
Concrete circular trajectory:
(1.00,0.00)→(0.98,0.17)→(0.76,0.64)→・・・・(radius*cosθ,radius*sinθ) 0≤θ≤360

Programmatically expressed as ....
for i in range(0,365,1):
 r = math.radians(i)
 x = radius1 * math.cos(r)
 y = radius1 * math.sin(r)
 print(x,y)
I think x,y will be circular.

I would also like to input the seismic displacement trajectory.
In the future, I would like to control this cylinder model so that it is equal to the input seismic displacement.
Specifically, I would like to control the cylinder model so that it becomes the same displacement by inputting the XY coordinate displacement for each earthquake time.
We would like to make the cylinder model oscillate in the XY direction by inputting the XY displacement waveform, the XY velocity waveform, or the XY displacement waveform of the earthquake.

To summarize
Is there any way/idea to move O.bodies.append(facet()) circularly?
Is there any idea how to make the cylinder shake in XY direction (like an earthquake)?

I apologize for my poor English.
Please ask questions if you do not understand.

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
内山康太郎
Solved:
Last query:
Last reply:
Revision history for this message
内山康太郎 (kenkoutaro) said : 		#1

#Current Program
from __future__ import print_function
from yade import pack,plot,polyhedra_utils,geom
from yade import export,qt
from yade.gridpfacet import *
import math
import numpy as np
import os

frictangle = np.radians(25)#45do
density = 3400.0
young = 3e8
gravity = (0.0, 0.0, 0)
velocity_cyl=5
velocity_cyl_idsCyl2=5
velocity_cyl_Box=5

# create cloud of spheres and insert them into the simulation
# we give corners, mean radius, radius variation

###############################################
#Create sphere
###############################################

mat_sp = CohFrictMat(alphaKr=.2,alphaKtw=.2,young = young, poisson = 0.35,frictionAngle=(frictangle), density=density)
O.materials.append(mat_sp)

pred=pack.inCylinder(centerBottom=(0,0,0.0),centerTop=(0,0,.02),radius=.03)
sp0=pack.randomDensePack(pred,spheresInCell=3000,radius=.0004)
O.bodies.append(sp0)

###############################################
#Boundary creation
###############################################

idsCyl1 = O.bodies.append(geom.facetCylinder((0, 0, 0.005), radius=.03, height=.01, segmentsNumber=30, wallMask=6))

idsCyl2 = O.bodies.append(geom.facetCylinder((0, 0, .015), radius=.03, height=.01, segmentsNumber=30, wallMask=5))

#idsCyl3 = O.bodies.append(geom.facetCylinder((0, 0, .0225), radius=0.10, height=.005, segmentsNumber=100, wallMask=2))

radius1 = .03
radius2 = .15

fixBoxIds=[]
moveBoxIds=[]

i=0

for i in range(0,365,5):
 r = math.radians(i)
 x1 = radius1 * math.cos(r)
 y1 = radius1 * math.sin(r)
 x2 = radius2 * math.cos(r-math.radians(5))
 y2 = radius2 * math.sin(r-math.radians(5))
 x3 = radius2 * math.cos(r+math.radians(5))
 y3 = radius2 * math.sin(r+math.radians(5))

 x4 = radius1 * math.cos(r-math.radians(5))
 y4 = radius1 * math.sin(r-math.radians(5))
 x5 = radius1 * math.cos(r+math.radians(5))
 y5 = radius1 * math.sin(r+math.radians(5))
 x6 = radius2 * math.cos(r)
 y6 = radius2 * math.sin(r)

 fixBoxIds.append(O.bodies.append(facet([(x1,y1,.01),(x2,y2,.01),(x3,y3,.01)])))
 fixBoxIds.append(O.bodies.append(facet([(x4,y4,.01),(x5,y5,.01),(x6,y6,.01)])))
 moveBoxIds.append(O.bodies.append(facet([(x1,y1,.01),(x2,y2,.01),(x3,y3,.01)])))
 moveBoxIds.append(O.bodies.append(facet([(x4,y4,.01),(x5,y5,.01),(x6,y6,.01)])))

upper_plate=box((0,0,-0.001),(0.06,0.06,0.001))
upper_plateID=O.bodies.append(upper_plate)
O.bodies[upper_plateID].state.blockedDOFs = 'xyXY'

#O.forces.addF(upper_plateID,(0,0,50e3),permanent=True)
O.forces.setPermF(upper_plateID,(0,0,100e3))
###############################################
#engines
###############################################

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                # interaction loop
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
                [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=gravity,damping=0.2),
        ##TranslationEngine(translationAxis=[0,0,1],velocity=-1,ids=idsCyl3),
        TranslationEngine(translationAxis=[1,0,0],velocity=velocity_cyl_idsCyl2,ids=idsCyl2),
        TranslationEngine(translationAxis=[1,0,0],velocity=velocity_cyl_Box,ids=moveBoxIds),
        PyRunner(command='checkStress()', iterPeriod=10),
        #PyRunner(command='checkStress_NO()', iterPeriod=10),
        PyRunner(command='plotPlotData()',iterPeriod=1),
]

def checkStress():
 cylDisplacement_x = O.bodies[idsCyl2[0]].state.displ()[0]
 cylDisplacement_y = O.bodies[idsCyl2[0]].state.displ()[1]
 plot.addData(
 cylDisplacement_x=cylDisplacement_x,
 cylDisplacement_y=cylDisplacement_y,
 )
 print('{:1.04f}'.format(cylDisplacement_x)," ", '{:1.04f}'.format(cylDisplacement_y))

###############################################
#Definition of distortion
###############################################
def checkStress_1():
 if O.bodies[idsCyl2[0]].state.displ()[0] > 0.020:
  O.pause()

def checkStress_2():
 if O.bodies[idsCyl2[0]].state.displ()[1] > 0.020:
  O.pause()

def checkStress_3():
 if O.bodies[idsCyl2[0]].state.displ()[0] < -0.020:
  O.pause()

def checkStress_4():
 if O.bodies[idsCyl2[0]].state.displ()[1] < -0.020:
  O.pause()

def checkStress_5():
 if O.bodies[idsCyl2[0]].state.displ()[0] > 0.0020:
  O.pause()

###############################################
#output
###############################################
sp_num = len([b for b in O.bodies if isinstance(b.shape, Sphere)]) # 粒子の総数をカウント
print("number of spheres = ",sp_num)

O.dt=0.5*PWaveTimeStep()
def plotPlotData():
 pipe_force_x = sum([O.forces.f(i)[0] for i in idsCyl2])
 pipe_force_y = sum([O.forces.f(i)[1] for i in idsCyl2])
 cylDisplacementx = O.bodies[idsCyl2[0]].state.displ()[0]
 cylDisplacementy = O.bodies[idsCyl2[0]].state.displ()[1]
 cylDisplacementz = O.bodies[idsCyl2[0]].state.displ()[2]
 Dz=upper_plate.state.displ()[2]
 plot.addData(
 i=O.iter,
 disp=O.time*velocity_cyl,
 force_x=pipe_force_x,
 force_y=pipe_force_y,
 Dz=Dz,
 cylDisplacementx=cylDisplacementx,
 cylDisplacementy=cylDisplacementy,
 cylDisplacementz=cylDisplacementz,
 )
 plot.saveDataTxt('pipe_force.out',vars=('disp','force_x','force_y',
      'Dz',
      'cylDisplacementx','cylDisplacementx','cylDisplacementx'
      ))

plot.plots={
 'disp':('force_x','force_y'),
 'i':('Dz'),
 'cylDisplacementx':('cylDisplacementy',)
 }
plot.plot()

#plot.saveDataTxt(O.tags['id'] + '.txt')

Revision history for this message
Jan Stránský (honzik) said : 		#2

Hello,

yes, you can use TranslationEngine to simulate any velocityVector(time) function, and therefore any trajectory.
Use PyRunner and set translationAxis and velocity of the TranslationEngine:

###
O.engines = [
    TranslationEngine(...,label="translation"),
    PyRunner(iterPeriod=1,command="setTranslation()"),
]
def setTranslation():
    t = O.time
    a = ... # any function of time
    v = ... # any function of time
    translation.translationAxis = a
    translation.velocity = v
###

Specifically for circular motion (with variable velocity):
###
r = 5
sph = sphere((r,0,0),1)
O.bodies.append(sph)
O.engines = [
    ForceResetter(),
    TranslationEngine(ids=[0],label="translation"),
    PyRunner(iterPeriod=1,command="setTranslation()"),
    NewtonIntegrator(),
]
O.dt = 1e-4

def setTranslation():
    x,y,z = sph.state.pos
    angle = atan2(y,x)
    v = sin(angle) + 1.4
    a = Vector3(-sin(angle),cos(angle),0)
    translation.translationAxis = a
    translation.velocity = v

view = yade.qt.View()
view.eyePosition = (0,0,20)
view.lookAt = (0,0,0)
###

Note that here a tangent of sphere is specified at each step, so technically the result is not exact circle, but "little" spiral.

Cheers
Jan

Revision history for this message
内山康太郎 (kenkoutaro) said : 		#3

Thank you very much Jan.

I want to move " idsCyl2 = O.bodies.append(geom.facetCylinder((0, 0, .015), radius=.03, height=.01, segmentsNumber=30, wallMask=5)) " with "radius 0.05" and "speed 1"

How can I do this?

Revision history for this message
Jan Stránský (honzik) said : 		#4

Just adapt the example in previous answer to your needs:
###
O.engines = [
    TranslationEngine(..., ids=idsCyl2, label="translation"),
    PyRunner(iterPeriod=1,command="setTranslation()"),
]
...
def setTranslation():
    ...
    v = 1
    ...
###
Cheers
Jan

Revision history for this message
内山康太郎 (kenkoutaro) said : 		#5

Thank you very much.
I post the completed program (excerpt).

##############################################################################
from __future__ import print_function
from yade import pack,plot,polyhedra_utils,geom
from yade import export,qt
from yade.gridpfacet import *
import math
import numpy as np
import os

idsCyl2 = O.bodies.append(geom.facetCylinder((0, 0, 0.015), radius=.06, height=.01, segmentsNumber=30, wallMask=5))

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                # interaction loop
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
                [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(),
        ##TranslationEngine(translationAxis=[0,0,1],velocity=-1,ids=idsCyl3),
        TranslationEngine(ids=idsCyl2,label="translation"),
        PyRunner(iterPeriod=1,command="setTranslation()"),
        PyRunner(command='checkStress()', iterPeriod=10),
        PyRunner(command='plotPlotData()',iterPeriod=1),
]

def setTranslation():
    x = O.bodies[idsCyl2[0]].state.displ()[0]
    y = O.bodies[idsCyl2[0]].state.displ()[1]
    angle = atan2(y,x)
    a = Vector3(-sin(angle),cos(angle),0)
    v = 0.001
    translation.translationAxis = a
    translation.velocity = v

def checkStress():
 cylDisplacement_x = O.bodies[idsCyl2[0]].state.displ()[0]
 cylDisplacement_y = O.bodies[idsCyl2[0]].state.displ()[1]
 plot.addData(
 cylDisplacement_x=cylDisplacement_x,
 cylDisplacement_y=cylDisplacement_y,
 )
 print('{:1.04f}'.format(cylDisplacement_x)," ", '{:1.04f}'.format(cylDisplacement_y))

O.dt=0.5*PWaveTimeStep()
def plotPlotData():
 cylDisplacementx = O.bodies[idsCyl2[0]].state.displ()[0]
 cylDisplacementy = O.bodies[idsCyl2[0]].state.displ()[1]
 plot.addData(
 i=O.iter,
 cylDisplacementx=cylDisplacementx,
 cylDisplacementy=cylDisplacementy,
 )
 plot.saveDataTxt('pipe_force.out',vars=('disp','force_x','force_y',
      'Dz',
      'cylDisplacementx','cylDisplacementx',
      ))

plot.plots={
 'disp':('force_x','force_y'),
 'i':('Dz'),
 'cylDisplacementx':('cylDisplacementy',)
 }
plot.plot()