New question #707913 on Yade: https://answers.launchpad.net/yade/+question/707913
Hi! In this simulation, the stress of the top plate(u1) is not increasing in the consolidation phase(before the shearing) although a velocity is imposed. Inside the shear box there is a compact parking of spheres. Can you please tell me if I am doing something wrong? # -*- coding: utf-8 -*- from yade import ymport, utils, plot,pack ################## Parâmetros ## copy spheres from the packing into the scene D_CILINDRO = 59.92e-3 H_CILINDRO = 26.8e-3 X_CUBO = D_CILINDRO Y_CUBO = H_CILINDRO Z_CUBO = D_CILINDRO # Parâmetros das partículas RAIO = 8e-4 VAR_RAIO = 10e-2 # 10% do Raio ESFERAS_CELULA = 2000 COR = (0,0,1) CORPO_DE_PROVA = 'cubo' PACKING = 'cilindro_0.0009_0.1_2000.sqlite' DAMPING = 0.5 dtCOEFF = 0.5 normalSTRESS = 50e3 # Pa normalVEL = 1e-3 shearVEL = 10e-5*X_CUBO # try different values to ensure quasi-static conditions intR=1#1.263 for X1Y1Z1_5k DENS=1600 YOUNG=40e6 FRICT=40 ALPHA=0.3 TENS=45e5 COH=45e6 iterMAX=40000 Output='shearTest_'+PACKING+'_K10_D3000_Y50e9A03F18T45e5C45e6_500kPa_shearVel1' ###################### Import of the sphere assembly ### material definition def sphereMat(): return CohFrictMat( young = YOUNG, frictionAngle = radians(FRICT), normalCohesion = 1, shearCohesion = 1, density = DENS, fragile = False) def wallMat(): return CohFrictMat( young = YOUNG, frictionAngle = 0, density = DENS, fragile = False) pred = pack.inAlignedBox((0,0,0),(X_CUBO,Y_CUBO, Z_CUBO)) #O.cell.hSize = Matrix3(D_CILINDRO, 0, 0, # 0,D_CILINDRO, 0, # 0, 0, H_CILINDRO) sp = SpherePack() memoizeDb = f'{CORPO_DE_PROVA}_{RAIO}_{VAR_RAIO}_{ESFERAS_CELULA}.sqlite' sp = pack.randomDensePack(pred, spheresInCell=2000, radius=RAIO, rRelFuzz = VAR_RAIO, memoizeDb = memoizeDb, returnSpherePack=True) spheres = sp.toSimulation(color=(0, 0, 1), material = sphereMat) ## preprocessing to get dimensions dim=utils.aabbExtrema() xinf=dim[0][0] xsup=dim[1][0] X=xsup-xinf yinf=dim[0][1] ysup=dim[1][1] Y=ysup-yinf zinf=dim[0][2] zsup=dim[1][2] Z=zsup-zinf ## initial surface S0=X*Z ## spheres factory R=0 Rmax=0 numSpheres=0. for o in O.bodies: if isinstance(o.shape,Sphere): o.shape.color=(0.7,0.5,0.3) numSpheres+=1 R+=o.shape.radius if o.shape.radius>Rmax: Rmax=o.shape.radius Rmean=R/numSpheres ###################### Criação da caixa cisalhante thickness=Y/100 oversizeFactor=1.3 ### Placas de carregamento ### loading platens O.bodies.append(utils.box(center=(xinf+X/2.,yinf-thickness+Rmean/10.,zinf+Z/2.), extents=(oversizeFactor*X/2,thickness,oversizeFactor*Z/2), material=wallMat,fixed=True)) l1 =O.bodies[-1] O.bodies.append(utils.box(center=(xinf+X/2.,ysup+thickness-Rmean/10.,zinf+Z/2.), extents=(oversizeFactor*X/2,thickness,oversizeFactor*Z/2), material=wallMat,fixed=True)) u1=O.bodies[-1] ### add top box O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,3*(ysup/4)+2*Rmean,zsup/2), extents=(thickness,Y/4,oversizeFactor*Z/2), material=wallMat, fixed=True)) u2=O.bodies[-1] O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,3*(ysup/4)+2*Rmean,zsup/2), extents=(thickness,Y/4,oversizeFactor*Z/2), material=wallMat, fixed=True)) u3=O.bodies[-1] O.bodies.append(utils.box(center=(xsup/2,3*(ysup/4)+2*Rmean,zinf-thickness+Rmean/10), extents=(oversizeFactor*X/2,Y/4,thickness), material=wallMat, fixed=True, wire=True)) u4=O.bodies[-1] O.bodies.append(utils.box(center=(xsup/2,3*(ysup/4)+2*Rmean,zsup+thickness-Rmean/10), extents=(oversizeFactor*X/2,Y/4,thickness), material=wallMat, fixed=True, wire=True)) u5=O.bodies[-1] ### add bottom box O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,(ysup/4)-2*Rmean,zsup/2), extents=(thickness,Y/4,oversizeFactor*Z/2), material=wallMat, fixed=True)) l2=O.bodies[-1] O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,(ysup/4)-2*Rmean,zsup/2), extents=(thickness,Y/4,oversizeFactor*Z/2), material=wallMat, fixed=True)) l3=O.bodies[-1] O.bodies.append(utils.box(center=(xsup/2,(ysup/4)-2*Rmean,zinf-thickness+Rmean/10), extents=(oversizeFactor*X/2,Y/4,thickness), material=wallMat, fixed=True, wire=True)) l4=O.bodies[-1] O.bodies.append(utils.box(center=(xsup/2,(ysup/4)-2*Rmean,zsup+thickness-Rmean/10), extents=(oversizeFactor*X/2,Y/4,thickness), material=wallMat, fixed=True, wire=True)) l5=O.bodies[-1] ###################### Engines O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='aabb'), Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=intR, label='ss2d3dg'), Ig2_Box_Sphere_ScGeom6D()], [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='interactionPhys')], [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(label='interactionLaw')] ), GlobalStiffnessTimeStepper(defaultDt=0.1*PWaveTimeStep(), timestepSafetyCoefficient=dtCOEFF), TranslationEngine(ids=[u1.id], translationAxis=(0.,-1.,0.), velocity=0., label='u1_yTranslation'), TranslationEngine(ids=[u1.id,u2.id,u3.id,u4.id,u5.id], translationAxis=(-1.,0.,0.), velocity=0., label='u_xTranslation'), TranslationEngine(ids=[l1.id,l2.id,l3.id,l4.id,l5.id], translationAxis=(1.,0.,0.), velocity=0., label='l_xTranslation'), PyRunner(iterPeriod=1, initRun=True, command='servoController()', label='servo'), NewtonIntegrator(damping=DAMPING, gravity=(0,-9.81,0), label='damper'), PyRunner(iterPeriod=int(iterMAX/1000), initRun=True, command='dataCollector()'), ] ###################### Engines definition ( servoController() and dataCollector() ) shearing=False sigmaN=0 tau=0 Fs1=0 Fs2=0 Xdispl=0 px0=0 Ydispl=0 py0=l1.state.pos[1] px0=l2.state.pos[0] py0=l1.state.pos[1] prevTranslation=0 n=0 F0 = normalSTRESS def servoController(): global px0, py0, sigmaN, n, Fn1, Fn2, shearing,u1, u2,u3,l1, l2, l3,Xdispl, Ydispl, tau,F0 sigmaN= abs(O.forces.f(u1.id)[1])/S0 # necessary? Xdispl = abs(l2.state.pos[0]-px0) tau = 0.0 placas = [u1,u2,u3,u4,u5,l1,l2,l3,l4,l5] for placa in placas: tau += abs(O.forces.f(placa.id)[0]) tau = tau/(2*S0) if abs(sigmaN-normalSTRESS)/normalSTRESS > 0.01: if sigmaN < normalSTRESS: u1_yTranslation.velocity += normalVEL/10 if sigmaN>(normalSTRESS): u1_yTranslation.velocity -= normalVEL/10 if abs((normalSTRESS-sigmaN)/normalSTRESS) < 0.01 and unbalancedForce() < 15/100: shearing = True print(f"sigmaN (Pa) = {sigmaN:.2f};", f"dif percentual (%) = {(normalSTRESS-sigmaN)/normalSTRESS*100.0:.2f};" f"velocidade da placa (m/s) = {u1.state.vel[1]:.3f};", f"unbalancedForce (%) = {unbalancedForce()*100.0:.1f};", f"tau (Pa) = {tau:.5f}") if shearing: u_xTranslation.velocity = shearVEL l_xTranslation.velocity = shearVEL def dataCollector(): global Xdispl, Ydispl, tau,sigmaN Xdispl = abs(l2.state.pos[0]-px0) Ydispl = abs(l1.state.pos[1]-py0) yade.plot.addData({'t':O.time, 'i':O.iter, 'Xdispl':Xdispl, 'sigmaN':sigmaN, 'tau':tau, 'Ydispl':Ydispl, 'unbF':utils.unbalancedForce()}) plot.saveDataTxt(Output) # defines window plot plot.plots={'i':('sigmaN','tau')} plot.plot() ################# graphical intervace from yade import qt qt.Controller() qt.View() v=qt.Renderer() v.dispScale=(1,1,1) # displacements are scaled for visualization purpose ################# to manage interaction detection factor during the first timestep O.step() ################# initializes the interaction detection factor to its default value (new contacts will be created between strictly contacting particles only) ss2d3dg.interactionDetectionFactor=-1. aabb.aabbEnlargeFactor=-1. ################# simulation starts here O.run(iterMAX) -- 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