Question #697421 on Yade changed: https://answers.launchpad.net/yade/+question/697421
Zhicheng Gao gave more information on the question: this is my code: ##______________ First section, generate sample_________ from __future__ import print_function from yade import pack, qt, plot from math import * nRead=readParamsFromTable( ## model parameters num_spheres=100, targetPorosity= .387, confiningPressure=-100000, ## material parameters compFricDegree=15,#contact friction during the confining phase finalFricDegree=30,#contact friction during the deviatoric loading young=2e8, poisson=.2, density=2600, alphaKr=7.5, alphaKtw=0, competaRoll=.22, finaletaRoll=.22, etaTwist=0, normalCohesion=0, shearCohesion=0, ## fluid parameters fluidDensity=1000, dynamicViscosity=.001, ## control parameters damp=0, stabilityThreshold=.001, ## output specifications filename='suffusion', unknowOk=True ) from yade.params.table import * O.periodic=True O.cell.hSize=Matrix3(.001,0,0, 0,.001,0, 0,0,.001) # create materials for spheres #shear strength is the sum of friction and adhesion, so the momentRotationLaw=True O.materials.append(CohFrictMat(alphaKr=alphaKr,alphaKtw=alphaKtw,density=density,etaRoll=competaRoll,etaTwist=etaTwist,frictionAngle=radians(compFricDegree),momentRotationLaw=True,normalCohesion=normalCohesion,poisson=poisson,shearCohesion=shearCohesion,young=young,label='spheres')) # generate particles packing sp=pack.SpherePack() sp.makeCloud((0,0,0),(.001,.001,.001),psdSizes=[0.00008,0.000125,0.0001592,0.0002003,0.0003153,0.000399,0.000502,0.0005743],psdCumm=[0.0,0.00628,0.0565,0.198,0.721,0.915,0.991,1.0],num=num_spheres,seed=1) sp.toSimulation(material='spheres') O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D()], [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='contact',setCohesionNow=False,setCohesionOnNewContacts=False)], [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True)], ), PeriodicFlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), PeriTriaxController(label='triax', # specify target values and whether they are strains or stresses goal=(confiningPressure,confiningPressure,confiningPressure), stressMask=7, # type of servo-control, the strain rate isn't determined, it shloud check the unbalanced force dynCell=True,maxStrainRate=(10,10,10), # wait until the unbalanced force goes below this value maxUnbalanced=stabilityThreshold,relStressTol=1e-3, doneHook='compactionFinished()' ), NewtonIntegrator(damping=0) ] qt.View() import sys def compactionFinished(): #check the current porosity # if the current porosity is lager than target Porosity and comFricDegree is lager than 10, # then we decrease friction value and apply it to all the bodies and contacts, # else we decrease rolling friction value. global compFricDegree, competaRoll if porosity()>targetPorosity and compFricDegree>5: # we decrease friction value and apply it to all the bodies and contacts compFricDegree = 0.95*compFricDegree setContactFriction(radians(compFricDegree)) print('Friction:', compFricDegree,'porosity:', porosity()) # python syntax, make each step printout sys.stdout.flush() elif porosity()>targetPorosity: # we decrease rolling fiction value and apply it to all the bodies and contacts competaRoll=0.95*competaRoll for b in O.bodies: b.mat.etaRoll=competaRoll for i in O.interactions: i.phys.etaRoll=competaRoll print('Rollingfriction:', b.mat.etaRoll, 'porosity:', porosity()) sys.stdout.flush() else: # after sample preparation, save the state O.save('compactedState'+filename+'.yade.gz') print('Compacted state saved', 'porosity', porosity()) # next time, called python command triax.doneHook='' O.pause() # enable energy tracking in the code O.trackEnergy=True # define function to record history def history(): plot.addData(unbalanced=unbalancedForce(),i=O.iter,exx=-triax.strain[0], eyy=-triax.strain[1], ezz=-triax.strain[2], sxx=-triax.stress[0],syy=-triax.stress[1],szz=-triax.stress[2], ev=-triax.strain[0]-triax.strain[1]-triax.strain[2], porosity=porosity(),Etot=O.energy.total(),**O.energy# save all available energy data ) O.engines=O.engines+[PyRunner(command='history()', iterPeriod=20)] # define what to plot plot.plots={'i':('unbalanced','porosity'),'i ':('sxx','syy','szz'),' i':('exx','eyy','ezz'),' i ':('Etot')} # show the plot plot.plot() plot.live=True O.run() O.wait() plot.saveDataTxt('compactedState'+filename) ##__________________________________second section, deviatoric loading__________________________ # change the contact parameters to the final calibration value setContactFriction(radians(finalFricDegree)) for b in O.bodies: b.mat.etaRoll=finaletaRoll for i in O.interactions: i.phys.etaRoll=finaletaRoll print(O.cell.hSize,O.cell.trsf) # set the current cell configuration to be the reference one O.cell.trsf=Matrix3.Identity print(O.cell.trsf, triax.strain) # change control type: keep constant confinement in x,y, 30% compression in z triax.goal=(confiningPressure,confiningPressure,-.3) triax.stressMask=3 # allow faster deformation along x,y to better maintain stresses triax.maxStrainRate=(1.,1.,.01) # call triaxFinished instead of compactionFinished triax.doneHook='triaxFinished()' def triaxFinished(): O.save('loadedState'+filename+'.yade.gz') print('Finished') O.pause() # Reset all plot data; keep plots and labels intact plot.resetData() O.run(1000,True) plot.saveDataTxt('loadedState'+filename) -- 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