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