Commit b19c09aa authored by Janek Kozicki's avatar Janek Kozicki Committed by Anton Gladky
Browse files

Small progress with #228, part of examples.

parent d81270ce
Pipeline #437664593 passed with stages
in 398 minutes and 34 seconds
......@@ -5,12 +5,12 @@
# Script to validate the lubrication component of the spring-dashpot law (Law2_ScGeom_ViscElPhys_Basic)
#
# Reproduce N times the simpleBallRebound script where two particles collides, varying the imposed effective restitution coefficient as a function of the stokes number
# Allows to reproduce curves of effective restitution coefficient as a function of the stokes number.
# Allows to reproduce curves of effective restitution coefficient as a function of the stokes number.
#
# All the quantities are made dimensionless
#
# THERE IS STILL A PROBLEM IN THE RESULT: AN UNEXPECTED CHANGE OF BEHAVIOR AT ST = 2e-2
#
#
#########################################################################################################################################################################
from yade import plot
import numpy as np
......@@ -19,79 +19,96 @@ import numpy as np
## Parameters of the configuration
#
stokesRange = np.logspace(0,3,200) #Stokes range explored for the collisions created. The length of stokes range correspond to the number of collisions that will be modeled.
stokesRange = np.logspace(
0, 3, 200
) #Stokes range explored for the collisions created. The length of stokes range correspond to the number of collisions that will be modeled.
## Input parameters
# Particles
normalRestit = 0.9 # Imposed dry restitution coefficient
kappa = 1e-6 # Stiffness dimensionless number, ratio between the particle stiffness and the granular pressure. Should be greater to 1e-4 to be in the rigid grain limit.
partMass = 1. # Dimensionless mass
diameter = 1. # Dimensionless diameters of the spheres
relativeVel = 1. # Dimensionless initial relative velocity
delta = 1e-4*diameter # Dimensionless initial relative position
roughnessRatio = 1e-2 # Roughness particle scale, dimensionless
normalRestit = 0.9 # Imposed dry restitution coefficient
kappa = 1e-6 # Stiffness dimensionless number, ratio between the particle stiffness and the granular pressure. Should be greater to 1e-4 to be in the rigid grain limit.
partMass = 1. # Dimensionless mass
diameter = 1. # Dimensionless diameters of the spheres
relativeVel = 1. # Dimensionless initial relative velocity
delta = 1e-4 * diameter # Dimensionless initial relative position
roughnessRatio = 1e-2 # Roughness particle scale, dimensionless
#Fluid
densityRatio = 1. # Particle to fluid density ratio
Stokes = 100. # Stokes number associated with the collision
densityRatio = 1. # Particle to fluid density ratio
Stokes = 100. # Stokes number associated with the collision
## Calculated parameters
# Particles
partDensity = partMass/(pi/6.*pow(diameter,3.)) # Particle density from its mass
O.dt = kappa*diameter/relativeVel # Time step in order to have a maximum overlap of kappa diameter
tContact = O.dt*300. # Contact time necessary to have a good resolution
partDensity = partMass / (pi / 6. * pow(diameter, 3.)) # Particle density from its mass
O.dt = kappa * diameter / relativeVel # Time step in order to have a maximum overlap of kappa diameter
tContact = O.dt * 300. # Contact time necessary to have a good resolution
# Fluid
fluidDensity = partDensity/densityRatio # Fluid density calculated to match the density ratio
fluidDensity = partDensity / densityRatio # Fluid density calculated to match the density ratio
gravityVector = Vector3(0.,0.,0.)
gravityVector = Vector3(0., 0., 0.)
#
## Creation of the framework
#
#Material
n = 0
effectiveRestit = np.zeros(len(stokesRange))
beadsCoupleID = np.zeros([len(stokesRange),2],dtype=int)
effectiveRestit = np.zeros(len(stokesRange))
beadsCoupleID = np.zeros([len(stokesRange), 2], dtype=int)
for Stokes in stokesRange:
nameMat = 'Mat'+str(n)
viscoDyn = partDensity*relativeVel*diameter/Stokes
effectiveRestit[n]=max(1e-3,normalRestit*(1. + 1./Stokes*log(2.*roughnessRatio)))
stiffnessContact = partMass/2./tContact*(log(effectiveRestit[n])**2+ pi**2) #See Schwager & Poschel 2007 on spring-dashpot contact law, estimated according contact stiffness
O.materials.append(ViscElMat(en=normalRestit, et = 1., frictionAngle = 0.5,young=stiffnessContact/diameter, poisson=0.5, density=partDensity,lubrication = True,viscoDyn = viscoDyn,roughnessScale=roughnessRatio*diameter, label=nameMat))
nameMat = 'Mat' + str(n)
viscoDyn = partDensity * relativeVel * diameter / Stokes
effectiveRestit[n] = max(1e-3, normalRestit * (1. + 1. / Stokes * log(2. * roughnessRatio)))
stiffnessContact = partMass / 2. / tContact * (
log(effectiveRestit[n])**2 + pi**2
) #See Schwager & Poschel 2007 on spring-dashpot contact law, estimated according contact stiffness
O.materials.append(
ViscElMat(
en=normalRestit,
et=1.,
frictionAngle=0.5,
young=stiffnessContact / diameter,
poisson=0.5,
density=partDensity,
lubrication=True,
viscoDyn=viscoDyn,
roughnessScale=roughnessRatio * diameter,
label=nameMat
)
)
#Create a free sphere at the origin with velocity in the direction of the other particle
beadsCoupleID[n,0] = O.bodies.append(sphere((2*n*diameter,0,0),diameter/2.0, material = nameMat))
beadsCoupleID[n,1] = O.bodies.append(sphere((2*n*diameter,0,diameter+ delta),diameter/2.0, material = nameMat))
for b in O.bodies:
if b.id ==beadsCoupleID[n,0]:
b.state.vel[2] = relativeVel/2.
if b.id ==beadsCoupleID[n,1]:
b.state.vel[2] = -relativeVel/2.
n+=1
beadsCoupleID[n, 0] = O.bodies.append(sphere((2 * n * diameter, 0, 0), diameter / 2.0, material=nameMat))
beadsCoupleID[n, 1] = O.bodies.append(sphere((2 * n * diameter, 0, diameter + delta), diameter / 2.0, material=nameMat))
for b in O.bodies:
if b.id == beadsCoupleID[n, 0]:
b.state.vel[2] = relativeVel / 2.
if b.id == beadsCoupleID[n, 1]:
b.state.vel[2] = -relativeVel / 2.
n += 1
#
## Simulation loop
#
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection'),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom() ],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
),
NewtonIntegrator(damping=0.0,gravity = gravityVector),
PyRunner(command = 'Meas()',iterPeriod = 100)
]
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(), Bo1_Facet_Aabb(), Bo1_Box_Aabb()], label='contactDetection'),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()], [Ip2_ViscElMat_ViscElMat_ViscElPhys()], [Law2_ScGeom_ViscElPhys_Basic()]
),
NewtonIntegrator(damping=0.0, gravity=gravityVector),
PyRunner(command='Meas()', iterPeriod=100)
]
O.saveTmp()
O.run()
effectiveRestitMeas = np.zeros(len(stokesRange))
stokesRangeMeas = np.zeros(len(stokesRange))
effectiveRestitMeas = np.zeros(len(stokesRange))
stokesRangeMeas = np.zeros(len(stokesRange))
def Meas():
if O.time > delta/relativeVel + tContact*500:
if O.time > delta / relativeVel + tContact * 500:
n = 0
k = 0
for b1 in O.bodies:
......@@ -99,18 +116,21 @@ def Meas():
n = 0
for idCouple in beadsCoupleID:
if b1.id == idCouple[0] and b2.id == idCouple[1]:
effectiveRestitMeas[k] = abs(b1.state.vel[2]-b2.state.vel[2])/relativeVel
effectiveRestitMeas[k] = abs(b1.state.vel[2] - b2.state.vel[2]) / relativeVel
stokesRangeMeas[k] = stokesRange[n]
k+=1
n+=1
k += 1
n += 1
print('finished, execute Plot() function now to see the results')
O.pause()
import matplotlib.pyplot as plt
def Plot():
plt.figure()
plt.semilogx(stokesRangeMeas,effectiveRestitMeas/normalRestit,'ok',label='measured')
plt.semilogx(stokesRange,effectiveRestit/normalRestit,'-r',label='imposed')
plt.semilogx(stokesRangeMeas, effectiveRestitMeas / normalRestit, 'ok', label='measured')
plt.semilogx(stokesRange, effectiveRestit / normalRestit, '-r', label='imposed')
plt.xlabel('Stokes number')
plt.ylabel('Effective restitution coefficient/dry restitution coefficient')
plt.show()
......@@ -2,14 +2,14 @@
# Author: Raphael Maurin, raphael.maurin@imft.fr
# 12/07/2021
#
# Script to validate the spring-dashpot law (Law2_ScGeom_ViscElPhys_Basic) and the lubrication force calculation from HydroForceEngine.
# Script to validate the spring-dashpot law (Law2_ScGeom_ViscElPhys_Basic) and the lubrication force calculation from HydroForceEngine.
#
# Consider a simple case of two spheres colliding with each other without gravity, each with a given opposite velocity of magnitude v/2
# Measure the ratio between the relative particles velocity after collision and the relative particles velocity before collision
# Measure the ratio between the relative particles velocity after collision and the relative particles velocity before collision
# Without lubrication or at Stokes > 10^3, this represents the dry restitution coefficient and one should recover the imposed value
#
# All the quantities are made dimensionless
#
#
#########################################################################################################################################################################
from yade import plot
import numpy as np
......@@ -22,71 +22,97 @@ lubricationForceApproach = False
## Input parameters
# Particles
normalRestit = 0.9 # Imposed dry restitution coefficient
kappa = 1e-6 # Stiffness dimensionless number, ratio between the particle stiffness and the granular pressure. Should be greater to 1e-4 to be in the rigid grain limit.
partMass = 1. # Dimensionless mass
diameter = 1. # Dimensionless diameters of the spheres
relativeVel = 1. # Dimensionless initial relative velocity
delta = 1e-3*diameter # Dimensionless initial relative position
roughnessRatio = 1e-2 # Roughness particle scale, dimensionless
normalRestit = 0.9 # Imposed dry restitution coefficient
kappa = 1e-6 # Stiffness dimensionless number, ratio between the particle stiffness and the granular pressure. Should be greater to 1e-4 to be in the rigid grain limit.
partMass = 1. # Dimensionless mass
diameter = 1. # Dimensionless diameters of the spheres
relativeVel = 1. # Dimensionless initial relative velocity
delta = 1e-3 * diameter # Dimensionless initial relative position
roughnessRatio = 1e-2 # Roughness particle scale, dimensionless
#Fluid
densityRatio = 100. # Particle to fluid density ratio
Stokes = 100. # Stokes number associated with the collision
effectiveRestit = normalRestit*(1. + 1./Stokes*log(2.*roughnessRatio))
densityRatio = 100. # Particle to fluid density ratio
Stokes = 100. # Stokes number associated with the collision
effectiveRestit = normalRestit * (1. + 1. / Stokes * log(2. * roughnessRatio))
## Calculated parameters
# Particles
partDensity = partMass/(pi/6.*pow(diameter,3.)) # Particle density from its mass
O.dt = kappa*diameter/relativeVel # Time step in order to have a maximum overlap of kappa diameter
tContact = O.dt*300. # Contact time necessary to have a good resolution
stiffnessContact = partMass/2./tContact*(log(effectiveRestit)**2+ pi**2) #See Schwager & Poschel 2007 on spring-dashpot contact law, estimated according contact stiffness
partDensity = partMass / (pi / 6. * pow(diameter, 3.)) # Particle density from its mass
O.dt = kappa * diameter / relativeVel # Time step in order to have a maximum overlap of kappa diameter
tContact = O.dt * 300. # Contact time necessary to have a good resolution
stiffnessContact = partMass / 2. / tContact * (
log(effectiveRestit)**2 + pi**2
) #See Schwager & Poschel 2007 on spring-dashpot contact law, estimated according contact stiffness
# Fluid
viscoDyn = partDensity*relativeVel*diameter/Stokes # Fluid viscosity calculated to match the stokes number required
fluidDensity = partDensity/densityRatio # Fluid density calculated to match the density ratio
viscoDyn = partDensity * relativeVel * diameter / Stokes # Fluid viscosity calculated to match the stokes number required
fluidDensity = partDensity / densityRatio # Fluid density calculated to match the density ratio
gravityVector = Vector3(0.,0.,0.)
gravityVector = Vector3(0., 0., 0.)
#
## Creation of the framework
#
#Material
if effectiveRestitutionApproach:
O.materials.append(ViscElMat(en=normalRestit, et = 1., frictionAngle = 0.4,young=stiffnessContact/diameter, poisson=0.5, density=partDensity,lubrication = True,viscoDyn = viscoDyn,roughnessScale=roughnessRatio*diameter, label='Mat1'))
O.materials.append(
ViscElMat(
en=normalRestit,
et=1.,
frictionAngle=0.4,
young=stiffnessContact / diameter,
poisson=0.5,
density=partDensity,
lubrication=True,
viscoDyn=viscoDyn,
roughnessScale=roughnessRatio * diameter,
label='Mat1'
)
)
else:
O.materials.append(ViscElMat(en=normalRestit, et = 1., frictionAngle = 0.4,young=stiffnessContact/diameter, poisson=0.5, density=partDensity, label='Mat1'))
O.materials.append(
ViscElMat(en=normalRestit, et=1., frictionAngle=0.4, young=stiffnessContact / diameter, poisson=0.5, density=partDensity, label='Mat1')
)
#Create a free sphere at the origin with velocity in the direction of the other particle
O.bodies.append(sphere((0,0,0),diameter/2.0, material = 'Mat1'))
O.bodies[0].state.vel[2] = relativeVel/2.
O.bodies.append(sphere((0, 0, 0), diameter / 2.0, material='Mat1'))
O.bodies[0].state.vel[2] = relativeVel / 2.
#Create a free sphere to send on the other sphere/the wall with a velocity v
O.bodies.append(sphere((0,0,diameter+ delta),diameter/2.0, material = 'Mat1'))
O.bodies[-1].state.vel[2] = -relativeVel/2.
O.bodies.append(sphere((0, 0, diameter + delta), diameter / 2.0, material='Mat1'))
O.bodies[-1].state.vel[2] = -relativeVel / 2.
#
## Simulation loop
#
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection'),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom() ],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
),
HydroForceEngine(densFluid = fluidDensity,viscoDyn = viscoDyn,zRef = 0.,gravity = gravityVector, deltaZ = diameter+ delta , nCell = 1,vCell = 1,ids = [0,1],label = 'hydroEngine',dead = True),
NewtonIntegrator(damping=0.0,gravity = gravityVector),
PyRunner(command = 'Plot()',iterPeriod = 100, label = 'plotFunction')
]
if lubricationForceApproach==True:
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(), Bo1_Facet_Aabb(), Bo1_Box_Aabb()], label='contactDetection'),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()], [Ip2_ViscElMat_ViscElMat_ViscElPhys()], [Law2_ScGeom_ViscElPhys_Basic()]
),
HydroForceEngine(
densFluid=fluidDensity,
viscoDyn=viscoDyn,
zRef=0.,
gravity=gravityVector,
deltaZ=diameter + delta,
nCell=1,
vCell=1,
ids=[0, 1],
label='hydroEngine',
dead=True
),
NewtonIntegrator(damping=0.0, gravity=gravityVector),
PyRunner(command='Plot()', iterPeriod=100, label='plotFunction')
]
if lubricationForceApproach == True:
hydroEngine.dead = False
hydroEngine.lubrication = True
hydroEngine.roughnessPartScale=roughnessRatio*diameter
hydroEngine.roughnessPartScale = roughnessRatio * diameter
hydroEngine.initialization()
O.saveTmp()
......@@ -94,13 +120,13 @@ O.run()
def Plot():
plot.addData(ratioVel = abs(O.bodies[1].state.vel[2]-O.bodies[0].state.vel[2])/relativeVel, dimensionlessTime = O.time/tContact)
if O.time > delta/relativeVel + tContact*500:
print("Effective restitution coefficient measured:",abs(O.bodies[1].state.vel[2]-O.bodies[0].state.vel[2])/relativeVel)
print("Effective restitution coefficient expected:",effectiveRestit)
plot.addData(ratioVel=abs(O.bodies[1].state.vel[2] - O.bodies[0].state.vel[2]) / relativeVel, dimensionlessTime=O.time / tContact)
if O.time > delta / relativeVel + tContact * 500:
print("Effective restitution coefficient measured:", abs(O.bodies[1].state.vel[2] - O.bodies[0].state.vel[2]) / relativeVel)
print("Effective restitution coefficient expected:", effectiveRestit)
O.pause()
#Plot the dimensionless sediment transport rate as a function of time during the simulation
plot.plots={'dimensionlessTime':('ratioVel')}
plot.plots = {'dimensionlessTime': ('ratioVel')}
plot.plot()
......@@ -6,60 +6,59 @@ Engine applying the linear motion of bicycle pedal e.g. moving points around the
from builtins import range
import time
## PhysicalParameters
Density=2400
frictionAngle=radians(35)
## PhysicalParameters
Density = 2400
frictionAngle = radians(35)
tc = 0.001
en = 0.3
es = 0.3
## Import wall's geometry
facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,tc=tc,en=en,et=es)) # **params sets kn, cn, ks, cs
sphereMat=O.materials.append(ViscElMat(density=Density,frictionAngle=frictionAngle,tc=tc,en=en,et=es))
facetMat = O.materials.append(ViscElMat(frictionAngle=frictionAngle, tc=tc, en=en, et=es)) # **params sets kn, cn, ks, cs
sphereMat = O.materials.append(ViscElMat(density=Density, frictionAngle=frictionAngle, tc=tc, en=en, et=es))
from yade import ymport
fctIds=O.bodies.append(ymport.stl('baraban.stl',color=(1,0,0),material=facetMat))
fctIds = O.bodies.append(ymport.stl('baraban.stl', color=(1, 0, 0), material=facetMat))
## Spheres
sphereRadius = 0.2
nbSpheres = (10,10,10)
nbSpheres = (10, 10, 10)
#nbSpheres = (50,50,50)
for i in range(nbSpheres[0]):
for j in range(nbSpheres[1]):
for k in range(nbSpheres[2]):
x = (i*2 - nbSpheres[0])*sphereRadius*1.1
y = (j*2 - nbSpheres[1])*sphereRadius*1.1
z = (k*2 - nbSpheres[2])*sphereRadius*1.1
s=sphere([x,y,z],sphereRadius,material=sphereMat)
O.bodies.append(s)
for j in range(nbSpheres[1]):
for k in range(nbSpheres[2]):
x = (i * 2 - nbSpheres[0]) * sphereRadius * 1.1
y = (j * 2 - nbSpheres[1]) * sphereRadius * 1.1
z = (k * 2 - nbSpheres[2]) * sphereRadius * 1.1
s = sphere([x, y, z], sphereRadius, material=sphereMat)
O.bodies.append(s)
## Timestep
O.dt=.2*tc
## Timestep
O.dt = .2 * tc
## Engines
O.engines=[
## Resets forces and momenta the act on bodies
ForceResetter(),
## Using bounding boxes find possible body collisions.
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
## Interactions
InteractionLoop(
## Create geometry information about each potential collision.
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
## Create physical information about the interaction.
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
## Constitutive law
[Law2_ScGeom_ViscElPhys_Basic()],
),
## Apply gravity
## Cundall damping must been disabled!
NewtonIntegrator(damping=0,gravity=[0,-9.81,0]),
## Saving results
#VTKRecorder(virtPeriod=0.04,fileName='/tmp/stlimp-',recorders=['spheres','facets']),
## Apply kinematics to walls
BicyclePedalEngine(ids=fctIds,rotationAxis=[0,0,1],radius = 1.2,angularVelocity=4.0)
## Engines
O.engines = [
## Resets forces and momenta the act on bodies
ForceResetter(),
## Using bounding boxes find possible body collisions.
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
## Interactions
InteractionLoop(
## Create geometry information about each potential collision.
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
## Create physical information about the interaction.
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
## Constitutive law
[Law2_ScGeom_ViscElPhys_Basic()],
),
## Apply gravity
## Cundall damping must been disabled!
NewtonIntegrator(damping=0, gravity=[0, -9.81, 0]),
## Saving results
#VTKRecorder(virtPeriod=0.04,fileName='/tmp/stlimp-',recorders=['spheres','facets']),
## Apply kinematics to walls
BicyclePedalEngine(ids=fctIds, rotationAxis=[0, 0, 1], radius=1.2, angularVelocity=4.0)
]
from yade import qt
qt.View()
#O.saveTmp()
#O.run()
......@@ -2,60 +2,59 @@
from builtins import range
import time
## PhysicalParameters
Density=2400
frictionAngle=radians(35)
## PhysicalParameters
Density = 2400
frictionAngle = radians(35)
tc = 0.001
en = 0.3
es = 0.3
## Import wall's geometry
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))
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))
from yade import ymport
fctIds=O.bodies.append(ymport.stl('baraban.stl',color=(1,0,0),material=facetMat))
fctIds = O.bodies.append(ymport.stl('baraban.stl', color=(1, 0, 0), material=facetMat))
## Spheres
sphereRadius = 0.2
nbSpheres = (10,10,10)
nbSpheres = (10, 10, 10)
#nbSpheres = (50,50,50)
for i in range(nbSpheres[0]):
for j in range(nbSpheres[1]):
for k in range(nbSpheres[2]):
x = (i*2 - nbSpheres[0])*sphereRadius*1.1
y = (j*2 - nbSpheres[1])*sphereRadius*1.1
z = (k*2 - nbSpheres[2])*sphereRadius*1.1
s=sphere([x,y,z],sphereRadius,material=sphereMat)
O.bodies.append(s)
for j in range(nbSpheres[1]):
for k in range(nbSpheres[2]):
x = (i * 2 - nbSpheres[0]) * sphereRadius * 1.1
y = (j * 2 - nbSpheres[1]) * sphereRadius * 1.1
z = (k * 2 - nbSpheres[2]) * sphereRadius * 1.1
s = sphere([x, y, z], sphereRadius, material=sphereMat)
O.bodies.append(s)
## Timestep
O.dt=.2*tc
## Timestep
O.dt = .2 * tc
## Engines
O.engines=[
## Resets forces and momenta the act on bodies
ForceResetter(),
## Using bounding boxes find possible body collisions.
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
## Interactions
InteractionLoop(
## Create geometry information about each potential collision.
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
## Create physical information about the interaction.
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
## Constitutive law
[Law2_ScGeom_ViscElPhys_Basic()],
),
## Apply gravity
## Cundall damping must been disabled!
NewtonIntegrator(damping=0,gravity=[0,-9.81,0]),
## Saving results
#VTKRecorder(virtPeriod=0.04,fileName='/tmp/stlimp-',recorders=['spheres','facets']),
## Apply kinematics to walls
RotationEngine(ids=fctIds,rotationAxis=[0,0,1],rotateAroundZero=True,angularVelocity=0.5)
## Engines
O.engines = [
## Resets forces and momenta the act on bodies
ForceResetter(),
## Using bounding boxes find possible body collisions.
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
## Interactions
InteractionLoop(
## Create geometry information about each potential collision.
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
## Create physical information about the interaction.
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
## Constitutive law
[Law2_ScGeom_ViscElPhys_Basic()],
),
## Apply gravity
## Cundall damping must been disabled!
NewtonIntegrator(damping=0, gravity=[0, -9.81, 0]),