------------------------------------------------------------ revno: 3728 committer: Christian Jakob <[email protected]> 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)
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
