New question #689117 on Yade: https://answers.launchpad.net/yade/+question/689117
Hello, I'm trying to obtain the permeability of spheres packed in a cylindrical shape using the PFV model in Yade. I previously posted a question about this problem [1]. The PVF only worked when I have my cylinder inside a box. The problem is that the box corners affects the flow results. Robert Caulk provided a solution to get the permeability in specific volume inside the cylinder [1, #20]. I tried that but the flow velocity results were high and unrealistic. I tried something that was suggested by Bruno Chareyre [1, #19] which is to block/deactivate the cells that are outside the cylinder core so that they're not included in the flow engine. I did that and I used Paraview to see the results and the cells around the cylinder were blocked, but the applied pressure seems to be only affecting the top surface [2, 3]. However, after the first ~1000 run, the cells were active again and the velocity suddenly jumps from 0.087 to 2269, you can see the Paraview visualization [4, 5]. Any idea why this happens and is there something wrong in my code (copied below)? Thank you, Othman P.S. to run the code, the sphere packing txt file can be found in [6]. [1] https://answers.launchpad.net/yade/+question/688685 [2] wireframe: https://drive.google.com/open?id=1nx5dLhYYd0g0hF-T-j74gs4mnXv9pSam [3] surface: https://drive.google.com/open?id=1Ftug3740yxHuuE0JJvsb71xafDoZJKUy [4] wireframe: https://drive.google.com/open?id=13e2Pq9pZE3phjqBgqKI7MswkfsMxMzit [5] surface: https://drive.google.com/open?id=1prgBJPlr5jJGrlhOnyj2q08Ot8s7TUzS [6] https://drive.google.com/open?id=17KEDSWCspG8SJtd4B3U32YxMhD-wzjm0 ------------------------------------------------ # -*- coding: utf-8 -*- from builtins import range from yade import pack, ymport import numpy as np import time tic=time.time() radiuscyl=.05 heightcyl=.25 dP=4.347e3 #Pa visc=1e-3 #Pa.sec - taken from Catalano, Chareyre, Bathelemy (2014) density=1000 #kg/m3 g=9.81 #m/s2 allx,ally,allz,r=np.genfromtxt('30Porcyl.txt', unpack=True) #access the packed cyl file (txt) mnx=min(allx)*0.99 #to get the xyz of the lowest corner, multiply by factor to reduce it to let the walls surround all the spheres. mny=min(ally)*0.99 mnz=min(allz)*0.99 mxx=max(allx)*1.01 #to get the xyz of the top corner, multiply by factor to increase it to let the wall surround all the spheres. mxy=max(ally)*1.01 mxz=max(allz)*1.01 print ('mn xyz, mx xyz', mnx,mny,mnz,mxx,mxy,mxz) 99998)) young=1e6 O.materials.append(FrictMat(young=young,poisson=0.2,density=1900,label='spheres')) O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls')) mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz) # corners of the box that contain the packign. Obtained from aabbExtrema walls=aabbWalls([mn,mx],thickness=0,material='walls') wallIds=O.bodies.append(walls) ############################# spheres ############################# spheres=O.bodies.append(ymport.text('/home/auser/DEM_codes/Flow/30Porcyl.txt')) Height=max(utils.aabbDim()) #Fix all particles in their positions. No deformation for i in O.bodies: i.state.blockedDOFs='xyzXYZ' O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop" ), FlowEngine(dead=1,label="flow"), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), NewtonIntegrator(damping=0.2) ] #Activate flow engine and set boundary conditions in order to get permeability flow.dead=0 flow.useSolver=3 #CHOLMOD: sparse direct solver for matrices. Its a library that is part of SuiteSparse linear algebra package. flow.permeabilityFactor=1 flow.viscosity=visc flow.bndCondIsPressure=[0,0,0,0,1,1] flow.bndCondValue=[0,0,0,0,0,dP] flow.boundaryUseMaxMin=[1,1,1,1,1,1] O.dt=1e-6 O.dynDt=False O.run(1,1) ################################################ blockcells ############################################################ ######################################################################################################################## numPoints = 100 xs = np.linspace(0.95*mnx,1.05*mxx,numPoints) ys = np.linspace(0.95*mny,1.05*mxy,numPoints) zs = np.linspace(0.95*mnz,1.05*mxz,numPoints) cellsHit = [] #get the point at the center of the cylinder cross-section cylcenterx=(mxx+mnx)/2 cylcentery=(mxy+mny)/2 for x,y,z in itertools.product(xs, ys, zs): cellId = flow.getCell(x,y,z) #getCell return cell id for a point in xyz if cellId in cellsHit: continue if np.sqrt((x-cylcenterx)**2+(y-cylcentery)**2) > radiuscyl: #this is to isolate the cells that are outside the cylinder. radiuscyl is the cylinder radius cellsHit.append(cellId) for i in cellsHit: # flow.blockCell(i,blockPressure=True) ##I tried this function but it didn't work and the cells in cellsHit were still counted in the flow model, so I used setCellPressure instead. flow.setCellPressure(i,0) O.run(1,1) qq=flow.averageVelocity() Q=abs(qq[2]) Permeability=abs((density*g*Q*Height)/dP) #hydraulic conductivity in distance/time print ("average velocity: ",Q, "Permeability m/s: ",Permeability,"Permeability in/hr: ", Permeability*141732) ######################################################################################################################## ######################################################################################################################## print ("runtime = ", time.time()-tic) #O.run() #O.engines=O.engines+[PyRunner(iterPeriod=20,command='flow.saveVtk(withBoundaries=True)')] -- You received this question notification because your team yade-users is an answer contact for Yade. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
