Question #703081 on Yade changed: https://answers.launchpad.net/yade/+question/703081
Status: Answered => Open Lukas Wasner is still having a problem: Hello Jan, thank you for the answer. Regarding to point [1], here is the "get to the right porosity" Script: readParamsFromTable(rMean=.0012, sig3=-100) from yade import pack, plot from yade.params.table import * mn,mx=Vector3(0,0,0),Vector3(.1,.1,.1) compFricDegree=30 damp=0.2 young=5e6 stabilityThreshold=0.01 targetPorosity=.4 O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2500,label='spheres')) O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls')) walls=aabbWalls([mn,mx],thickness=0,material='walls') wallIds=O.bodies.append(walls) sp=pack.SpherePack() sp.makeCloud(mn,mx,rMean=rMean,seed=1) sp.makeCloud(mn,mx,rMean=.5*rMean,seed=1) sp.toSimulation(material='spheres') triax=TriaxialStressController( thickness = 0, stressMask = 7, internalCompaction=False, ) newton=NewtonIntegrator(damping=damp) 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()] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), triax, #TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key), newton ] triax.goal1=triax.goal2=triax.goal3=sig3 while 1: O.run(1000, True) unb=unbalancedForce() print('unbalanced force:',unb,' mean stress: ',triax.meanStress) if unb<stabilityThreshold and abs(sig3-triax.meanStress)/sig3<0.001: print('Compression equilibrium reached') break import sys while triax.porosity>targetPorosity: compFricDegree*=.95 setContactFriction(radians(compFricDegree)) print('\r Friction: ',compFricDegree,' porosity:',triax.porosity,) sys.stdout.flush() O.run(250,1) if triax.porosity<targetPorosity: sp.fromSimulation() sp.save(f'packtest_rMean{rMean}_sig{sig3}_n_{triax.porosity}.txt') print('Pack ist exportiert') O.pause() waitIfBatch() --------------------------------------------------------------------- Regarding point [2]: I have tried the volume argument and unfortunately get different results: ------------------------------------------------------------------------------- For script 3 (Predicate is a box). I get a similar porosity when I insert the volume of the predicate into the porosity function. .5672 without argument .56898 with Predicate Volume as argument from yade import pack sp=SpherePack() sp.load('packtest_rMean0.0012_sig-100.txt') print('target porosity is 0.38542') pred=pack.inAlignedBox(minAABB=(-.01,-.01,-.01), maxAABB=(.01, .01, .01)) O.bodies.append(geom.facetBox((0, 0, 0), (.01, .01, .01), wallMask=31)) spf=pack.filterSpherePack(predicate=pred,spherePack=sp,returnSpherePack=True) spf.toSimulation() print('porosity without argument',porosity()) print('porosity with argument volume',porosity(volume=(.02*.02*.02))) --------------------------------------------------------------------------------------------------------------- For script 2 (Predicate is a cylinder). I get a lower porosity if I insert the volume of the predicate into the porosity function. Nevertheless, this is much too high. ---------------------------------------------------------------------------------------------------------------- from yade import pack sp=SpherePack() sp.load('packtest_rMean0.0012_sig-100.txt') print('target porosity is 0.38542') O.bodies.append(geom.facetCylinder(center=(0, 0, .01), radius=.01, height=.02,segmentsNumber=50, wallMask=6)) pred=pack.inCylinder(centerBottom=(0,0,0), centerTop=(0,0,.02), radius=.01) spf2=pack.filterSpherePack(predicate=pred,spherePack=sp,returnSpherePack=True) spf2.toSimulation() print('porosity without argument',porosity()) print('porosity with argument', porosity(volume=((pi*.01**2)*.02))) ------------------------------------------------------------------------------------------------------------ Unfortunately, at this point I do not know how to proceed, since now only the proof has been given that the porosity function uses the volume of the predicate, which it is supposed to do. Why there is a larger deviation with the cylinder predicate, I cannot explain myself. I would still like to find a way to work with the porosity function, since I would like to map the pore number over the stress later in the oedometer experiment. I conclude that the filtered section of the imported sphere packing has a smaller total particle volume in relation to the new (predicate) volume than the total particle volume in relation to the final compressed particle cloud (old volume) from [1]. Since in [1] the particle cloud was compressed isotropically, I assumed that the porosity is equally distributed over the entire volume and the cropped part accordingly also has the same porosity as the entire sphere packing volume. How do I manage to get the same porosity in my predicate as it is in the old volume (imported pack)? Thanks in advance, Lukas -- 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 : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp