Question #688791 on Yade changed: https://answers.launchpad.net/yade/+question/688791
Leonard posted a new comment: Hi Jan, Thanks for your suggestion. Here are the two scripts based on Yade 2018.02b: ######### The following is the script1 which is for generating sample (takes around 30 seconds): # -*- coding: utf-8 -*- #************************************************************************* from __future__ import division from yade import pack, plot import math import numpy as np import timeit ########## this script is modified from https://gitlab.com/yade- dev/trunk/blob/master/examples/triax-tutorial/script-session1.py#L142 utils.readParamsFromTable(lowerR=10.0,upperR=10.0) from yade.params import table num_spheres=2000# number of spheres targetPorosity = 0.5 #the porosity we want for the packing compFricDegree = 30 # initial contact friction during the confining phase (will be decreased during the REFD compaction process) finalFricDegree = 30 # contact friction during the deviatoric loading rate=-0.01 # loading rate (strain rate) damp=0.4 # damping coefficient stabilityThreshold=0.001 # we test unbalancedForce against this value in different loops (see below) young=5e7 # contact stiffness confinement=100e3 mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07) ## create materials for spheres and plates MatWall=O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls')) MatSand = O.materials.append(CohFrictMat(isCohesive=False,young=young,poisson=0.5,frictionAngle=radians(30), density=2650.0,normalCohesion=1e6, shearCohesion=1e6,label='sand')) ## create walls around the packing walls=aabbWalls([mn,mx],thickness=0,material='walls') wallIds=O.bodies.append(walls) ## use a SpherePack object to generate a random loose particles packing sp=pack.SpherePack() sp.makeCloud(mn,mx,-1,0.3333,num_spheres,True, 0.95,seed=1) O.bodies.append([sphere(center,rad,material='sand') for center,rad in sp]) Gl1_Sphere.quality=3 triax=TriaxialStressController( maxMultiplier=1.+6e5/young, # spheres growing factor (fast growth) finalMaxMultiplier=1.+2e4/young, # spheres growing factor (slow growth) thickness = 0, stressMask = 7, internalCompaction=True, # If true the confining pressure is generated by growing particles ) newton=NewtonIntegrator(damping=damp) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()], [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(), Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), triax, TriaxialStateRecorder(iterPeriod=100,file='WallStresses'), newton, ] #Display spheres with 2 colors for seeing rotations better Gl1_Sphere.stripes=0 #the value of (isotropic) confining stress defines the target stress to be applied in all three directions triax.goal1=triax.goal2=triax.goal3=-confinement while 1: O.run(1000, True) unb=unbalancedForce() print 'unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1] if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001: break print "### Isotropic completed ###" import sys while triax.porosity>targetPorosity: # we decrease friction value and apply it to all the bodies and contacts compFricDegree = 0.95*compFricDegree setContactFriction(radians(compFricDegree)) print "\r Friction: ",compFricDegree," porosity:",triax.porosity, sys.stdout.flush() O.run(500,1) print "### state 2 Reach target porosity completed ###" print 'top', -triax.stress(triax.wall_top_id)[1], 'right',-triax.stress(triax.wall_right_id)[0],'front',-triax.stress(triax.wall_front_id)[2] triax.internalCompaction=False triax.stressMask = 7 triax.goal1=triax.goal2=triax.goal3=-confinement triax.max_vel=0.001 O.save('sample.yade.gz') ##################################################### The following is script2 which is for triaxial test, it illustrates the problem. ## # -*- coding: utf-8 -*- #************************************************************************* from yade import pack, plot import math import numpy as np O.load('sample.yade.gz') Gl1_Sphere.quality=3 confinement=100e3 triax=TriaxialStressController( ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions. ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002 thickness = 0, stressMask = 5, internalCompaction=False, # If true the confining pressure is generated by growing particles ) newton=NewtonIntegrator(damping=0.4) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()], [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(), Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), triax, TriaxialStateRecorder(iterPeriod=100,file='WallStresses'), PyRunner(command='stopIfDamaged()',iterPeriod=1000), newton, ] Gl1_Sphere.stripes=0 yade.qt.Controller() yade.qt.View() setContactFriction(radians(30)) triax.stressMask = 5 triax.goal2=-0.015 triax.goal1=-confinement triax.goal3=-confinement from yade import plot ### a function saving variables def history(): plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2], Volumetric_strian=triax.strain[0]+triax.strain[1]+triax.strain[2], s11=-triax.stress(triax.wall_right_id)[0], s22=-triax.stress(triax.wall_top_id)[1], s33=-triax.stress(triax.wall_front_id)[2], dev=-triax.stress(triax.wall_top_id)[1]-confinement, i=O.iter) def stopIfDamaged(): if -triax.strain[1]>0.15: O.pause() plot.saveDataTxt('data') if 1: O.engines=O.engines[0:5]+[PyRunner(iterPeriod=1000,command='history()',label='recorder')]+[PyRunner(iterPeriod=100,command='stopIfDamaged()',label='checkDamage')]+O.engines[5:7] plot.plots={'e22':('s11','s22','s33'),'e22 ':('dev'),'e22 ':('Volumetric_strian')} plot.labels={'e22':'$\epsilon_{a}$','s11':'$\sigma_{11}$','s22':'$\sigma_{22}$','s33':'$\sigma_{33}$'} plot.plot() ######## Run script2 after script1, I have: In [1]: FATAL /build/yade-fDuCoe/yade-2018.02b/core/ThreadRunner.cpp:30 run: Exception occured: InsertionSortCollider.verletDist>0, but unable to locate NewtonIntegrator within O.engines. You may switch the position of newton with PyRunner in the O.engines (in script2) to see the difference. For me, the error is solved by this way. Cheers, Leonard -- 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
