New question #687838 on Yade: https://answers.launchpad.net/yade/+question/687838
Hi everyone, I am using Ubuntu 18.04, and Yade 2019-08-08.git-775ae74 I use the Triaxial code by Bruno Chareyre [1] to understand how it works. [1] https://gitlab.com/yade-dev/trunk/blob/master/examples/triax-tutorial/script-session1.py My goal is to code a Triaxial test with 20000 particles in different sizes (sand) under different strain and stress paths and get all inter-particle forces, branch vectors, contact normals, and fabric tensor. (with the properties defined in the code) This is the code I am running. My questions are at the end of this message. ################################################### ############################################################################################################################ ######### TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z IS THE FRONT AXIS ######### ############################################################################################################################ import datetime aa = datetime.datetime.now() print ('************** START **************') import numpy as np import math from yade import pack, plot, qt, export, utils from datetime import datetime ###################################################### ######### DEFINING VARIABLES ######### print ('============ DEFINING VARIABLES ============') nRead=readParamsFromTable( num_spheres=20000, compFricDegree = 30, key='_triax_base_', unknownOk=True ) from yade.params import table num_spheres=table.num_spheres key=table.key targetPorosity = 0.35 compFricDegree = table.compFricDegree finalFricDegree = 30 rate=-0.02 damp=0.2 stabilityThreshold=0.01 young=5e6 poisson=0.3 mn,mx=Vector3(0,0,0),Vector3(20,20,20) sigmaIso=-50e3 ###################################################### ######### DEFINING MATERIALS ######### print ('============ DEFINING MATERIALS ============') O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=2600,label='spheres')) O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='frictionlesswalls')) #################################################### ######### DEFINING PACKING ######### print ('============ DEFINING PACKING ============') frictionlesswalls=aabbWalls([mn,mx],thickness=0,material='frictionlesswalls') wallIds=O.bodies.append(frictionlesswalls) sp=pack.SpherePack() clumps=False if clumps: volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2]) mean_rad = pow(0.09*volume/num_spheres,0.3333) c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)]) sp.makeClumpCloud(mn,mx,[c1],periodic=False) sp.toSimulation(material='spheres') O.bodies.updateClumpProperties() else: volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2]) sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp]) ########################################################## ######### DEFINING TRIAXIAL TEST ######### print ('============ DEFINING TRIAXIAL TEST ============') triax=TriaxialStressController( maxMultiplier=1.+2e4/young, finalMaxMultiplier=1.+2e3/young, thickness = 0, stressMask = 7, internalCompaction=True, ) newton=NewtonIntegrator(damping=damp) #################################################### ######### DEFINING ENGINES ######### print ('============ DEFINING ENGINES ============') 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()] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), triax, TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key), newton ] Gl1_Sphere.stripes=True if nRead==0: yade.qt.Controller(), yade.qt.View() print ('Number of elements: ', len(O.bodies)) print ('Box Volume: ', triax.boxVolume) print ('Box Volume calculated: ', volume) ############################################################### ######### APPLYING CONFINING PRESSURE ######### print ('============ APPLYING CONFINING PRESSURE ============') triax.stressmask=7 triax.goal1=sigmaIso triax.goal2=sigmaIso triax.goal3=sigmaIso while 1: O.run(1000, True) unb=unbalancedForce() meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3 print ('unbalanced force:',unb,' mean stress engine: ',triax.meanStress,' mean stress (Calculated): ',meanS) print ('porosity=',triax.porosity) print ('void ratio=',triax.porosity/(1-triax.porosity)) print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~') if unb<stabilityThreshold and abs(sigmaIso-triax.meanStress)/sigmaIso<0.001: break O.save('confinedPhase'+key+'.xml') print ('################## Isotropic phase is finished and saved successfully ##################') e22Check=triax.strain[1] ###################################################### ######### DEVIATORIC LOADING ######### print ('============ APPLYING DEVIATORIC LOADING ============') triax.internalCompaction=False setContactFriction(radians(finalFricDegree)) triax.stressMask = 5 triax.goal1=sigmaIso triax.goal2=rate triax.goal3=sigmaIso newton.damping=0.1 unb=unbalancedForce() axialS=triax.stress(triax.wall_top_id)[1] print ('step=', O.iter, 'unbalanced force:',unb,' sigma2: ',axialS, 'q=', axialS-sigmaIso) print ('axial deformation (%)', (triax.strain[1]-e22Check)*100) print ('~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~') O.save('final.xml') print ('################## Deviatoric phase is finished and saved successfully ##################') ###################################################### ######### DEFINING FUNCTIONS ######### print ('============ DEFINING FUNCTIONS ============') def history(): plot.addData(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], i=O.iter) ######################################################## ######### RECORD AND PLOT DATA ######### print ('============ RECORD AND PLOT DATA ============') O.run(5000,True) plot.plots={'e22':('s11','s22','s33',None,'ev')} plot.labels={'s11':'$\sigma_{11}$' , 's22':'$\sigma_{22}$' , 's33':'$\sigma_{33}$' , 'e22':'$\epsilon_{22}$' , 'ev':'$\epsilon_{V}$'} plot.plot() plot.saveDataTxt('results'+key) plot.saveGnuplot('plotScript'+key) print ('************** END **************') bb = datetime.datetime.now() cc = bb - aa print( int(cc.total_seconds())) print(elapsed_time) ################################################### Question: 1- When the simulation is finished? When I run the code, It takes some minutes to show the ************** END ************** message on the screen, while the Play button has not been chosed in the Yade window yet. If I press the Play button, it takes some hours and then I get run errors. 2- In the code we have defined the strain rate for the deviatoric part, what's the criterion to finish the test? does it work with damping? I have read about it but it's not clear for me yet. 3- I want to get the time of running but I get the error for it. I have used datetime function for it as it is in the code. 4- I don't get the conventional graph of triaxial test. the graph window shows up but it is just a point. 5- If I get the simple triaxial test correctly, how can I get the fabric tensor, inter-particle forces, branch vectors, and contact normals as the raw data? Thank you. -- 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 : [email protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
