yade-users team mailing list archive

[Jan Stránský <question660155@xxxxxxxxxxxxxxxxxxxxx>]

Question #660155 on Yade changed:
https://answers.launchpad.net/yade/+question/660155

Jan Stránský proposed the following answer:
Hi Rodolfo,

based on personal communication, please find attached an example.
Drag coefficient can be also implemented as tensor if needed.
Before relying on this engine, please properly test it.(!!!)

You have to define area yourself, in local coordinate system (with respect to principal axes of inertia). But it also can be automated.
If performance would be a problem, C++ implementation is also possible.

cheers
Jan

####
from yade import plot

class PyDragEngine:
  def __init__(self,rho,cd,area,id):
    """
    Python implementation of 
    rho ... mass density of fluid
    cd ... drag coefficient
    area ... "area tensor" in local coordinate system
    id ... body id to apply force on

    TODO cd could also be a tensor (as area)
    TODO testing
    """
    self.rho = float(rho)
    self.cd = float(cd)
    self.area = Matrix3(area)
    self.id = int(id)
  def __call__(self):
    b = O.bodies[self.id]
    v = b.state.vel.norm() # velocity magnitude
    if abs(v) < 1e-6: # do nothing for very low velocity
      return
    n = b.state.vel / v # velocity direction as unit vector
    r = b.state.ori.toRotationMatrix() # rotation matrix
    a = r*self.area*r.transpose() # rotation of area tensor according to rot of body. Maybe it is rT*a*r, please check
    a = n.dot(a*n) # projection of "a" to direction of n.
    f = .5*self.rho*self.cd*a*pow(v,2) # drag force magnitude
    f = -f*n # drag force vector
    O.forces.addF(self.id,f) # add it to the body

cid,sids = O.bodies.appendClumped(( # testing clump, ax=2, ay=3, az=6
  sphere((0,0,0),.5),
  sphere((1,0,0),.5),
  sphere((2,0,0),.5),
  sphere((0,1,0),.5),
  sphere((1,1,0),.5),
  sphere((2,1,0),.5),
))
clump = O.bodies[cid]
O.dt = 0
clump.state.ori = Quaternion((1,0,0),.5*pi)*Quaternion((0,0,1),.5*pi) # to test with non-unit ori, then ax=3, ay=6, az=2

a = Matrix3(3,0,0, 0,6,0, 0,0,2) # global area tensor, diag(ax,ay,az)
r = clump.state.ori.toRotationMatrix()
a = r.transpose()*a*r # local area tensor # ?? r*a*rT ??
pyDragEngine = PyDragEngine(rho=1.225,cd=.6,area=a,id=cid)

O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb()]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom()],
    [Ip2_FrictMat_FrictMat_FrictPhys()],
    [Law2_ScGeom_FrictPhys_CundallStrack()]
  ),
  PyRunner(iterPeriod=1,command="pyDragEngine()"), # before NowtonIntegrator (!)
  NewtonIntegrator(damping=0.,label='newton'),
  PyRunner(iterPeriod=1,command="plotAddData()"),
]

def plotAddData():
  plot.addData(t=O.time,v=clump.state.vel.norm())

# test for different gravity directions
O.dt = 3e0
newton.gravity = Vector3(1,0,0) # ax=3
O.run(100,True)
calm()
newton.gravity = Vector3(0,1,0) # ay=6
O.run(100,True)
calm()
newton.gravity = Vector3(0,0,1) # az=2
O.run(100,True)
calm()
newton.gravity = Vector3(1,1,1).normalized()
O.run(100,True)

plot.plots = {'t':'v'}
plot.plot()
####

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.