yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10149
[Branch ~yade-pkg/yade/git-trunk] Rev 3721: fix buoyancy example for clumps
------------------------------------------------------------
revno: 3721
committer: Christian Jakob <jakob@xxxxxxxxxxxxxxxxxxx>
timestamp: Thu 2013-10-17 11:11:53 +0200
message:
fix buoyancy example for clumps
modified:
examples/clumps/apply-buoyancy-clumps.py
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'examples/clumps/apply-buoyancy-clumps.py'
--- examples/clumps/apply-buoyancy-clumps.py 2013-10-04 08:11:01 +0000
+++ examples/clumps/apply-buoyancy-clumps.py 2013-10-17 09:11:53 +0000
@@ -18,9 +18,9 @@
F_buo = -volumeOfDisplacedWater*fluidDensity*gravityAcceleration.'''
#define material properties:
-shear_modulus = 3.2e10
-poisson_ratio = 0.15
-young_modulus = 2*shear_modulus*(1+poisson_ratio)
+shearModulus = 3.2e10
+poissonRatio = 0.15
+youngModulus = 2*shearModulus*(1+poissonRatio)
angle = 0.5 #friction angle in radians
rho_p = 2650 #density of particles
rho_f = 5000 #density of the fluid rho_f > rho_p = floating particles
@@ -38,7 +38,7 @@
waterColor = (.2,.2,.7)#blue
#material:
-id_Mat=O.materials.append(FrictMat(young=young_modulus,poisson=poisson_ratio,density=rho_p,frictionAngle=angle))
+id_Mat=O.materials.append(FrictMat(young=youngModulus,poisson=poissonRatio,density=rho_p,frictionAngle=angle))
Mat=O.materials[id_Mat]
#create assembly of spheres:
@@ -55,7 +55,7 @@
b.shape.color=clumpColor
#create boundary:
-O.bodies.append(facetBox((0,0,1), (boundaryMax-boundaryMin)/2, fixed=True, material=Mat, color=boxColor))#boundaryMax-boundaryMin
+O.bodies.append(geom.facetBox((0,0,1), (boundaryMax-boundaryMin)/2, fixed=True, material=Mat, color=boxColor))
#define engines:
O.engines=[
@@ -63,68 +63,54 @@
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
- [Ip2_FrictMat_FrictMat_MindlinPhys(betan=0.0,betas=0.0,label='ContactModel')],
- [Law2_ScGeom_MindlinPhys_Mindlin(neverErase=False,label='Mindlin')]
+ [Ip2_FrictMat_FrictMat_MindlinPhys()],
+ [Law2_ScGeom_MindlinPhys_Mindlin(neverErase=False)]
),
- GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8, defaultDt=0.9*PWaveTimeStep(),label='ts'),
NewtonIntegrator(gravity=(0,0,-10),damping=0.7,label='integrator')
]
-ts.dead=True
#definition to apply buoyancy:
-def apply_buo(water_height,saturatedList):
+def applyBuoyancy():
+ global waterLevel
+ waterLevel = (O.time-t0) * v_iw + boundaryMin[2]
for b in O.bodies:
- if b not in saturatedList:
- h_low = 1e9
- h_high = -1e9
- if b.isStandalone and isinstance(b.shape,Sphere):
- rad = b.shape.radius
- pos = b.state.pos
- h_low = pos[2] - rad
- h_high = pos[2] + rad
- elif b.isClump: #determine rad, h_high and h_low for clumps:
- keys = O.bodies[b.id].shape.members.keys()
- for ii in range(0,len(keys)):
- pos = O.bodies[keys[ii]].state.pos
- h_low = min(h_low,pos[2]-O.bodies[keys[ii]].shape.radius)
- h_high = max(h_high,pos[2]+O.bodies[keys[ii]].shape.radius)
- rad = ( 3*b.state.mass/(4*pi*O.bodies[keys[0]].mat.density) )**(1./3.) #get equivalent radius from clump mass
- else:
- continue
- if water_height > h_low:
- dh = water_height - h_low
- dh = min(dh,2*rad) #to get sure, that dh is not bigger than 2*radius
- F_buo = -1*(pi/3)*dh*dh*(3*rad - dh)*rho_f*integrator.gravity # = -V*rho*g
- #apply buoyancy force (will overwrite old forces applied with addF command)
- if b.isStandalone and isinstance(b.shape,Sphere):
- O.forces.addF(b.id,F_buo,permanent=True)
- if b.isClump:
- keys = O.bodies[b.id].shape.members.keys()
- for ii in range(0,len(keys)):
- O.forces.addF(keys[ii],(O.bodies[keys[ii]].state.mass/b.state.mass)*F_buo,permanent=True) # F_buo_clumpmember = massPortion*F_buo_clump
- if water_height > h_high:
- saturatedList.append(b)
+ zMin = 1e9
+ zMax = -1e9
+ if b.isStandalone and isinstance(b.shape,Sphere):
+ rad = b.shape.radius
+ zMin = b.state.pos[2] - rad
+ dh = min((waterLevel - zMin),2*rad) #to get sure, that dh is not bigger than 2*radius
+ elif b.isClump: #determine rad, zMin and zMax for clumps:
+ for ii in b.shape.members.keys():
+ pos = O.bodies[ii].state.pos
+ zMin = min(zMin,pos[2]-O.bodies[ii].shape.radius)
+ zMax = max(zMax,pos[2]+O.bodies[ii].shape.radius)
+ #get equivalent radius from clump mass:
+ rad = ( 3*b.state.mass/(4*pi*O.bodies[b.shape.members.keys()[0]].mat.density) )**(1./3.)
+ #get dh relative to equivalent sphere, but acting when waterLevel is between clumps z-dimensions zMin and zMax:
+ dh = min((waterLevel - zMin)*2*rad/(zMax - zMin),2*rad)
else:
- if water_height < h_high:
- saturatedList.remove(b) #remove "hoppers" from saturatedList
- return saturatedList
+ continue
+ if dh > 0:
+ F_buo = -1*(pi/3)*dh*dh*(3*rad - dh)*rho_f*integrator.gravity # = -V*rho*g
+ O.forces.addF(b.id,F_buo,permanent=True)
#STEP1: reduce overlaps from replaceByClumps() method:
O.dt=1e-6 #small time step for preparation steps via calm()
print '\nSTEP1 in progress. Please wait a minute ...\n'
O.engines=O.engines+[PyRunner(iterPeriod=10000,command='calm()',label='calmRunner')]
-O.run(1000000,True)
+O.run(100000,True)
+
+#STEP2: let particles settle down
calmRunner.dead=True
-
-#STEP2: let particles settle down
-O.dt=1e-5
+O.dt=2e-5
print '\nSTEP2 in progress. Please wait a minute ...\n'
O.run(50000,True)
#start PyRunner engine to apply buoyancy:
-saturatedList = []
t0 = O.time
-O.engines=O.engines+[PyRunner(iterPeriod=100,command='saturatedList=apply_buo(((O.time-t0) * v_iw),saturatedList)',label='buoLabel')]
+waterLevel = boundaryMin[2]
+O.engines=O.engines+[PyRunner(iterPeriod=100,command='applyBuoyancy()',label='buoLabel')]
#create 2 facets, that show water height:
idsWaterFacets = []
@@ -146,5 +132,5 @@
v=qt.View()
v.eyePosition=(-7,0,2); v.upVector=(0,0,1); v.viewDir=(1,0,-.1); v.axes=True; v.sceneRadius=1.9
print '\nSTEP3 started ...\n'
+
O.run(70000)
-