yade-dev team mailing list archive

[ChiaWeng Boon <chiaweng@xxxxxxxxx>]

Hello everyone,

If you run the attached script long enough, you will see kinetic
energy increasing.  There is no damping so it should be constant.   Do
you get what I get too? Do you think somewhere went wrong?


Regards,

Boon

# booncw
# chiaweng@xxxxxxxxx 

#!/usr/local/bin/yade-trunk -x
# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-

# Script to test a sphere-sphere interaction in displacement control

from yade import utils
import math
import random

# encoding: utf-8
from yade import pack

sp=pack.SpherePack()
## corners of the initial packing
mn,mx=Vector3(-10,-10,-10),Vector3(10,10,10)

# makeCloud parameters "documented" only by the argument names in the c++ signature now:
# http://beta.arcig.cz/~eudoxos/yade/epydoc/yade._packSpheres.SpherePack-class.html#makeCloud
## box between mn and mx, avg radius .5 ± ½(.5*.2), 10k spheres (will be less, obviously), not periodic
sp.makeCloud(mn,mx,.2,0.0,3000,False)

## create material #0, which will be used as default
O.materials.append(FrictMat(young=150e6,poisson=.4,frictionAngle=0.0,density=2600))
O.materials.append(FrictMat(young=150e6,poisson=.4,frictionAngle=0.0,density=2600,label='frictionless'))


## copy spheres from the packing into the scene
## use default material, don't care about that for now
O.bodies.append([utils.sphere(s[0],s[1],material='frictionless') for s in sp])
for b in O.bodies:
	v=[random.random(),random.random(),random.random()]
	normV = 1.0*sqrt(v[0]**2+v[1]**2+v[2]**2)
	signV1 = random.random()	
	signV2 = random.random()
	signV3 = random.random()
	if signV1 < 0.5:
		signV1 =-1.0
	else:
		signV1 = 1.0
	if signV2< 0.5:
		signV2 =-1.0
	else:
		signV2 = 1.0
	if signV3 < 0.5:
		signV3 =-1.0
	else:
		signV3 = 1.0
	b.state.vel=[signV1*v[0]/normV,signV2*v[1]/normV,signV3*v[2]/normV]
	b.state.vel=5.0*b.state.vel

#__________________________________________________________________
# list of engines

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],nBins=5,sweepLength=.05),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack(label='c')]
	),
	NewtonIntegrator(damping=0.0,exactAsphericalRot=False),
	PyRunner(iterPeriod=1000,command='myAddPlotData()')
]



#__________________________________________________________________

#__________________________________________________________________
from yade import qt
qt.View()
qt.Controller()

#__________________________________________________________________
# plot some results
from math import *
from yade import plot







O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away


def myAddPlotData():
	uf=utils.unbalancedForce()
	KE2=utils.kineticEnergy()
	plot.addData(t1=O.iter,t2=O.iter,t3=O.iter,kineticEnUtils=KE2)


plot.plots={'t1':('kineticEnUtils')}


O.bodies.append(utils.box(center=(0,0,-10),extents=(9.9,9.9,0.05),orientation=[1,0,0,0],wire=True,dynamic=False,material='frictionless'))
O.bodies.append(utils.box(center=(0,0,10),extents=(9.9,9.9,0.05),orientation=[1,0,0,0],wire=True,dynamic=False,material='frictionless'))
O.bodies.append(utils.box(center=(10,0,0),extents=(0.05,9.9,9.9),orientation=[1,0,0,0],wire=True,dynamic=False,material='frictionless'))
O.bodies.append(utils.box(center=(-10,0,0),extents=(0.05,9.9,9.9),orientation=[1,0,0,0],dynamic=False,material='frictionless',wire=True))
O.bodies.append(utils.box(center=(0,10,0),extents=(9.9,0.05,9.9),orientation=[1,0,0,0],dynamic=False,material='frictionless',wire=True))
O.bodies.append(utils.box(center=(0,-10,0),extents=(9.9,0.05,9.9),orientation=[1,0,0,0],dynamic=False,material='frictionless',wire=True))

Attachment: KEwithTime.png
Description: PNG image