Question #669048 on Yade changed: https://answers.launchpad.net/yade/+question/669048
SayedHessam posted a new comment: Dear Robert, Sorry for the lack of information. BTW, you can find herewith my script, as you requested: ############################ ### DEFINING ENGINES ### ############################ 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 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth) finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth) thickness = 0, ## switch stress/strain control using a bitmask. What is a bitmask, huh?! ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0. ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4 ## to put it differently, the mask is the integer whose binary representation is xyz, i.e. ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc. 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_ScGeom(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()] ), ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf) FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), triax, TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key), newton ] #Display spheres with 2 colors for seeing rotations better Gl1_Sphere.stripes=0 if nRead==0: yade.qt.Controller(), yade.qt.View() ## UNCOMMENT THE FOLLOWING SECTIONS ONE BY ONE ## DEPENDING ON YOUR EDITOR, IT COULD BE DONE ## BY SELECTING THE CODE BLOCKS BETWEEN THE SUBTITLES ## AND PRESSING CTRL+SHIFT+D ####################################### ### APPLYING CONFINING PRESSURE ### ####################################### #the value of (isotropic) confining stress defines the target stress to be applied in all three directions triax.goal1=triax.goal2=triax.goal3=-10000 while 1: O.run(1000, True) #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium unb=unbalancedForce() print 'unbalanced force:',unb,' mean stress: ',triax.meanStress if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001: break O.save('confinedState'+key+'.yade.gz') print "### Isotropic state saved ###" print triax.porosity print triax.meanStress print len(O.bodies) ##################################################### ### Example of how to record and plot data ### ##################################################### from yade import plot ## a function saving variables def history(): plot.addData(unbalanced=unbalancedForce(),e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2], ev=-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], devi = -triax.stress(triax.wall_top_id)[1] - (-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2]) / 2.0, p = triax.meanStress, i=O.iter) if 1: # include a periodic engine calling that function in the simulation loop O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7] #O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder')) else: # With the line above, we are recording some variables twice. We could in fact replace the previous # TriaxialRecorder # by our periodic engine. Uncomment the following line: O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder') O.run(100,True) ## declare what is to plot. "None" is for separating y and y2 axis plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')} ## the traditional triaxial curves would be more like this: #plot.plots={'e22':('s11','s22','s33',None,'ev')} # display on the screen (doesn't work on VMware image it seems) plot.plot() devi_test = -triax.stress(triax.wall_top_id)[1] - (-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2]) / 2.0 def stress_control1(): global devi_test while devi_test < 50000: O.run(100) devi_test = -triax.stress(triax.wall_top_id)[1] - (-triax.stress(triax.wall_right_id)[0]\ -triax.stress(triax.wall_front_id)[2]) / 2.0 def stress_control2(): global devi_test while devi_test > -50000: O.run(100) devi_test = -triax.stress(triax.wall_top_id)[1] - (-triax.stress(triax.wall_right_id)[0]\ -triax.stress(triax.wall_front_id)[2]) / 2.0 #B. Activate flow engine and set boundary conditions in order to get permeability flow.dead=0 flow.defTolerance=0.3 flow.meshUpdateInterval=200 flow.useSolver=3 flow.permeabilityFactor=1 flow.viscosity=1.3e-3 flow.fluidBulkModulus = 2e9 flow.imposePressure(Vector3(triax.width/2,triax.height/2,triax.depth/2),0.001) flow.bndCondIsPressure=[0,0,0,0,0,0] flow.bndCondValue=[0,0,0,0,0,0] flow.boundaryUseMaxMin=[0,0,0,0,0,0] O.dt=0.1e-3 O.dynDt=False O.run(1,1) Qin = flow.getBoundaryFlux(2) Qout = flow.getBoundaryFlux(3) tic_toc = 0 print(O.iter) print(triax.stressMask) print(triax.goal1) print(triax.goal2) print(triax.goal3) Regards Sam -- 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