ringSimpleViscoelastic.py 2.32 KB
Newer Older
Anton Gladky's avatar
Anton Gladky committed
1 2
#!/usr/bin/python
# -*- coding: utf-8 -*-
3

4
from __future__ import print_function
5
from builtins import range
6
from yade import ymport
7 8 9 10 11 12 13 14 15 16 17 18 19

## Omega
o=Omega() 

## PhysicalParameters 
Density=2400
frictionAngle=radians(35)
sphereRadius=0.05
tc = 0.001
en = 0.3
es = 0.3

## Import wall's geometry
20 21
facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,tc=tc,en=en,et=es))
sphereMat=O.materials.append(ViscElMat(density=Density,frictionAngle=frictionAngle,tc=tc,en=en,et=es))
22 23

walls = O.bodies.append(ymport.stl('ring.stl',material=facetMat))
24 25 26

def fill_cylinder_with_spheres(sphereRadius,cylinderRadius,cylinderHeight,cylinderOrigin,cylinderSlope):
	spheresCount=0
27 28
	for h in range(0,int(cylinderHeight/sphereRadius/2)):
			for r in range(1,int(cylinderRadius/sphereRadius/2)):
29
				dfi = asin(0.5/r)*2
30
				for a in range(0,int(6.28/dfi)):
31 32 33
					x = cylinderOrigin[0]+2*r*sphereRadius*cos(dfi*a)
					y = cylinderOrigin[1]+2*r*sphereRadius*sin(dfi*a)
					z = cylinderOrigin[2]+h*2*sphereRadius
34
					s=sphere([x,y*cos(cylinderSlope)+z*sin(cylinderSlope),z*cos(cylinderSlope)-y*sin(cylinderSlope)],sphereRadius,material=sphereMat)
35 36 37 38 39 40 41
					o.bodies.append(s)
					spheresCount+=1
	return spheresCount

# Spheres
spheresCount=0
spheresCount+=fill_cylinder_with_spheres(sphereRadius,0.5,0.10,[0,0,0],radians(0))
42
print("Number of spheres: %d" % spheresCount)
43 44 45 46


## Engines 
o.engines=[
47
	## Resets forces and momenta the act on bodies
48
	ForceResetter(),
49

50 51
	## Using bounding boxes find possible body collisions.
	InsertionSortCollider([
Václav Šmilauer's avatar
Václav Šmilauer committed
52 53
		Bo1_Sphere_Aabb(),
		Bo1_Facet_Aabb(),
54
	]),
55
	# Interactions
56
	InteractionLoop(
57 58 59
		## Create geometry information about each potential collision.
		[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
		## Create physical information about the interaction.
60
		[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
61
		## Constitutive law
Anton Gladky's avatar
Anton Gladky committed
62
		[Law2_ScGeom_ViscElPhys_Basic()],
63
	),
64
	## Apply gravity
65
	## Cundall damping must been disabled!
66
	NewtonIntegrator(damping=0,gravity=[0,-9.81,0]),
67
	## Apply kinematics to walls
68
    ## angularVelocity = 0.73 rad/sec = 7 rpm
69
	RotationEngine(ids=walls,rotationAxis=[0,0,1],rotateAroundZero=True,angularVelocity=0.73)
70

71 72
]

73 74
for b in O.bodies:
	if isinstance(b.shape,Sphere):
75
		b.state.blockedDOFs='zXY'
76

77
O.dt=0.02*tc
78

79
O.saveTmp('init');
80 81 82

from yade import qt
renderer=qt.Renderer()
Anton Gladky's avatar
Anton Gladky committed
83
renderer.wire=True
84 85 86
#qt.Controller()
qt.View()
O.run()
87