New question #662105 on Yade: https://answers.launchpad.net/yade/+question/662105
Hello there, I have been struggling with this for a while now and could not figure out what is wrong. My simulation (shear test with periodic boundary conditions) explodes as soon as the flow engine is activated (see the last 2 lines) and I don't know if it is due to the technical problems that I am facing with the flow engine or if there is something wrong with my script. Could someone please try it and see if it works on his side? Many thanks Luc ### from yade import pack sp=pack.SpherePack() O.periodic=True # dimensions of sample RADIUS=0.05 length=10*(2*RADIUS) height=length width=length thickness=RADIUS/10. # microproperties DENS=2600 E=1e8 P=0.5 compFRIC=10. FRIC=30. TENS=0. COH=0. # loading conditions SN=5e6 RATE=0.1 #### O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width) O.materials.append(CohFrictMat(isCohesive=True,density=DENS,young=E,poisson=P,frictionAngle=radians(0.),normalCohesion=1e100,shearCohesion=1e100,label='boxMat')) O.materials.append(CohFrictMat(isCohesive=True,density=DENS,young=E,poisson=P,frictionAngle=radians(compFRIC),normalCohesion=TENS,shearCohesion=COH,label='sphereMat')) upBox = utils.box(center=(0,2*height+thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat') lowBox = utils.box(center=(0,height-thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat') O.bodies.append([upBox,lowBox]) sp.makeCloud((0,height+1.5*RADIUS,0),(length,2*height-1.5*RADIUS,width),rMean=RADIUS,rRelFuzz=0.2,periodic=True) O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='sphereMat') for s in sp]) effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2] volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol #### engines flow=PeriodicFlowEngine( isActivated=0 ,useSolver=3 ### what should we define in the next 2 lines? ,defTolerance=-1 # with this value, triangulation is done according to meshUpdateInterval ,meshUpdateInterval=1000 # this value by default ,duplicateThreshold=0.5 # how do we define this? ,boundaryUseMaxMin=[0,0,0,0,0,0] # what is this for? should we have [0,0,1,1,0,0] because of the top and bottom walls? ,wallIds=[-1,-1,1,0,-1,-1] # top wall id=0, bottom wall id=0 ,wallThickness=thickness ,bndCondIsPressure=[1,1,0,0,1,1] # drained in the horizontal plane ,bndCondValue=[0,0,0,0,0,0] ,permeabilityFactor=1.e-3 ,viscosity=1. ,label='flowEng' ) O.engines=[ ForceResetter() ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True) ,InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()], [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()], [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()] ) ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep()) ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-1.e5/volRatio,-1.e5/volRatio,-1.e5/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='compactionDone()',label='triax') ,flow ,NewtonIntegrator(damping=0.2,label='newton') ] def compactionDone(): print 'Changing contact properties now' tt=TriaxialCompressionEngine() tt.setContactProperties(FRIC) print 'DONE! iter=',O.iter,'| isotropic confinement done: stresses=',volRatio*triax.stress, triax.dead=True O.pause() #### Initialization print 'SAMPLE PREPARATION!' O.run(1000000,1) #### Applying normal stress print 'NORMAL LOADING! iter=',O.iter stage=0 stiff=fnPlaten=currentSN=0. def servo(): global stage,stiff,fnPlaten,currentSN if stage==0: currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]) unbF=unbalancedForce() boundaryVel=copysign(min(0.1,abs(0.5*(currentSN-SN))),currentSN-SN) O.bodies[0].state.vel[1]=boundaryVel if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ): stage+=1 fnPlaten=O.forces.f(0)[1] print 'Normal stress =',currentSN,' | unbF=',unbF for i in O.interactions.withBody(O.bodies[0].id): stiff+=i.phys.kn print 'DONE! iter=',O.iter O.pause() if stage==1: fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2]) boundaryVel=copysign(min(RATE,abs(0.35*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired) # why 0.35? O.bodies[0].state.vel[1]=boundaryVel O.engines = O.engines[:4]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[4:] O.run(1000000,1) print 'Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]) print 'Normal stress (contacts) = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1] #### preparing for shearing print 'Gluing spheres to boundary walls' ## gluing particles in contact with the walls for i in O.interactions: if i.isReal: if isinstance(O.bodies[i.id1].shape,Box): O.bodies[i.id2].shape.color[0]=1 if isinstance(O.bodies[i.id2].shape,Box): O.bodies[i.id1].shape.color[0]=1 for i in O.interactions: if i.isReal and ( O.bodies[i.id1].shape.color[0]==1 and O.bodies[i.id2].shape.color[0]==1 ): O.bodies[i.id1].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion O.bodies[i.id2].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion O.bodies[i.id1].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion O.bodies[i.id2].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion i.phys.initCohesion=True #### shearing print 'SHEARING! iter=',O.iter flow.isActivated=1 # sample explodes as soon as fluid is activated??? O.bodies[0].state.vel[0]=RATE O.run(10) -- 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