yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
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)
-