yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3728
committer: Christian Jakob <jakob@xxxxxxxxxxxxxxxxxxx>
timestamp: Fri 2013-10-04 10:11:01 +0200
message:
  adapt buoyancy example with new O.forces.addF() method and make it look nicer
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-06-24 13:23:34 +0000
+++ examples/clumps/apply-buoyancy-clumps.py	2013-10-04 08:11:01 +0000
@@ -4,19 +4,18 @@
 '''	The script shows how to include the effect of buoyancy in a particle assembly 
 	under condition of an increasing water-level. Water-level at start of the sim-
 	ulation is at the lower boundary of the model. During the calculation particles
-	gets partly surrounded by water and buoyancy (=density*volumeOfDisplacedWater)
+	get partly surrounded by water and buoyancy (=fluidDensity*volumeOfDisplacedWater)
 	is increasing until particle is completely surrounded. When particle is sur-
 	rounded volumeOfDisplacedWater is equal to the volume of the particle.
 	
-	For calculation of buoyancy of a clump the theoretical radius 
-	R = (3*mass/(4*pi*density))^1/3 
+	For calculation of buoyancy of a clump the equivalent radius 
+	R = (3*clumpMass/(4*pi*particleDensity))^1/3 
 	of a sphere with clump mass and clump volume 
-	V = (4*pi*R^3)/3 = mass/density
+	V = (4*pi*R^3)/3 = clumpMass/particleDensity
 	is used. This approximation can be used for well rounded clumps.
 	
-	Buoyancy is included via reduced mass (=massAtStart - dm), where dm is the 
-	buoyancy. Reduced mass must be corrected via density scaling for calculation
-	of correct particle accelerations.'''
+	Buoyancy is included with an additional force 
+	F_buo = -volumeOfDisplacedWater*fluidDensity*gravityAcceleration.'''
 
 #define material properties:
 shear_modulus			= 3.2e10
@@ -24,6 +23,7 @@
 young_modulus			= 2*shear_modulus*(1+poisson_ratio)
 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
 
 v_iw					= 1  #velocity of increasing water-level
 
@@ -31,6 +31,12 @@
 boundaryMin = Vector3(-1.5,-1.5,0)
 boundaryMax = Vector3(1.5,1.5,2)
 
+#define colors:
+sphereColor = (.8,.8,0.)#dirty yellow
+clumpColor	= (1.,.55,0.)#dark orange
+boxColor	= (.1,.5,.1)#green
+waterColor	= (.2,.2,.7)#blue
+
 #material:
 id_Mat=O.materials.append(FrictMat(young=young_modulus,poisson=poisson_ratio,density=rho_p,frictionAngle=angle))
 Mat=O.materials[id_Mat]
@@ -38,14 +44,18 @@
 #create assembly of spheres:
 sp=pack.SpherePack()
 sp.makeCloud(minCorner=boundaryMin,maxCorner=boundaryMax,rMean=.2,rRelFuzz=.5,num=100,periodic=False)
-O.bodies.append([sphere(c,r,material=Mat) for c,r in sp])
+O.bodies.append([sphere(c,r,material=Mat,color=sphereColor) for c,r in sp])
 
 #create template for clumps and replace 10% of spheres by clumps:
 templateList = [clumpTemplate(relRadii=[1,.5],relPositions=[[0,0,0],[.7,0,0]])]
 O.bodies.replaceByClumps(templateList,[.1])
+#color clumps:
+for b in O.bodies:
+	if b.isClumpMember:
+		b.shape.color=clumpColor
 
 #create boundary:
-O.bodies.append(facetBox((0,0,1), (boundaryMax-boundaryMin)/2, fixed=True, material=Mat))#boundaryMax-boundaryMin
+O.bodies.append(facetBox((0,0,1), (boundaryMax-boundaryMin)/2, fixed=True, material=Mat, color=boxColor))#boundaryMax-boundaryMin
 
 #define engines:
 O.engines=[
@@ -62,7 +72,7 @@
 ts.dead=True
 
 #definition to apply buoyancy:
-def apply_buo(water_height,saturatedList,startMass):
+def apply_buo(water_height,saturatedList):
 	for b in O.bodies:
 		if b not in saturatedList:
 			h_low = 1e9
@@ -78,18 +88,25 @@
 					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*startMass[b.id]/(4*pi*O.bodies[keys[0]].mat.density) )**(1./3.)		#get radius from startMass
+				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
-				dm = (pi/3)*dh*dh*(3*rad - dh)*1000		#buoyancy acts via reduced mass dm=rho*V_cap
-				b.state.mass = startMass[b.id] - dm		#needs a list of masses of all particles at start (see 5-final.py)
-				#reduced mass must be corrected via density scaling for accelerations a = F/m_reduced, a_correct = a*densityScaling, see line 179 in pkg/dem/NewtonIntegrator.cpp:
-				b.state.densityScaling = (startMass[b.id] - dm)/startMass[b.id]
+				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)
+		else:
+			if water_height < h_high:
+				saturatedList.remove(b)	#remove "hoppers" from saturatedList
 	return saturatedList
 
 #STEP1: reduce overlaps from replaceByClumps() method:
@@ -99,47 +116,35 @@
 O.run(1000000,True)
 calmRunner.dead=True
 
-#STEP2: let particles settle down:
+#STEP2: let particles settle down
 O.dt=1e-5
 print '\nSTEP2 in progress. Please wait a minute ...\n'
 O.run(50000,True)
 
-#get maximum body id:
-idmax=0
-for b in O.bodies:
-	idmax=max(idmax,b.id)
-integrator.densityScaling=True	#activate density scaling
-
-#get a list of masses at start:
-startMass = [0 for ii in range(0,idmax+1)]
-for b in O.bodies:
-	if (b.isStandalone and isinstance(b.shape,Sphere)) or b.isClump:
-		startMass[b.id] = b.state.mass
-
 #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,startMass)',label='buoLabel')]
+O.engines=O.engines+[PyRunner(iterPeriod=100,command='saturatedList=apply_buo(((O.time-t0) * v_iw),saturatedList)',label='buoLabel')]
 
 #create 2 facets, that show water height:
 idsWaterFacets =  []
 idsWaterFacets.append(O.bodies.append(facet( \
 			[ boundaryMin, [boundaryMax[0],boundaryMin[1],boundaryMin[2]], [boundaryMax[0],boundaryMax[1],boundaryMin[2]] ], \
-			fixed=True,material=FrictMat(young=0),color=(0,0,1),wire=False)))#no interactions will appear
+			fixed=True,material=FrictMat(young=0),color=waterColor,wire=False)))#no interactions will appear
 idsWaterFacets.append(O.bodies.append(facet( \
 			[ [boundaryMax[0],boundaryMax[1],boundaryMin[2]], [boundaryMin[0],boundaryMax[1],boundaryMin[2]], boundaryMin ], \
-			fixed=True,material=FrictMat(young=0),color=(0,0,1),wire=False)))#no interactions will appear
+			fixed=True,material=FrictMat(young=0),color=waterColor,wire=False)))#no interactions will appear
 
 #set velocity of incr. water level to the facets:
 for ids in idsWaterFacets:
 	O.bodies[ids].state.vel = [0,0,v_iw]
 
-#STEP3: simulate buoyancy#
+#STEP3: simulate buoyancy with increasing water table condition
 O.dt=3e-5
 from yade import qt
 qt.Controller()
 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(60000)
+O.run(70000)