New question #687871 on Yade: https://answers.launchpad.net/yade/+question/687871
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) Based on the code and its output which can be seen in the following, I have asked my questions at the end of this message. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ############################################################################################################################ ######### TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z IS THE FRONT AXIS ######### ############################################################################################################################ print ('************** START **************') import numpy as np import datetime import time import math from yade import 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 FUNCTIONS ######### print ('============ DEFINING FUNCTIONS ============') from yade import plot 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, n = triax.porosity) ###################################################### ######### 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) from yade import pack 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, ) #################################################### ######### 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, PyRunner(iterPeriod=20,command='history()',label='recorder'), TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key), NewtonIntegrator(damping=damp,label="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 ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~') 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)) 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 ##################') ######################################################## ######### 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 **************') O.realtime print ('Analysis has been taken for',O.realtime, 'seconds or', O.realtime/60, 'minutes') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Now output: ehsan@ehsan:~/Desktop$ /home/ehsan/yade/install/bin/yade-2019-08-08.git-775ae74 Q.py Welcome to Yade 2019-08-08.git-775ae74 Using python version: 3.6.9 (default, Nov 7 2019, 10:44:02) [GCC 8.3.0] TCP python prompt on localhost:9000, auth cookie `scdyua' XMLRPC info provider on http://localhost:21000 Running script Q.py ************** START ************** ============ DEFINING VARIABLES ============ ============ DEFINING FUNCTIONS ============ ============ DEFINING MATERIALS ============ ============ DEFINING PACKING ============ ============ DEFINING TRIAXIAL TEST ============ ============ DEFINING ENGINES ============ The constructor with a shareWidget is deprecated, use the regular contructor instead. Number of elements: 20006 Box Volume: 6.90392657805617e-310 Box Volume calculated: 8000.0 ============ APPLYING CONFINING PRESSURE ============ ~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~ unbalanced force: 0.007284146090760937 mean stress engine: -1.3899805976256234 mean stress (Calculated): -2.537923714531715 porosity= 0.8343868539188172 void ratio= 5.038168005756043 ################## Isotropic phase is finished and saved successfully ################## ============ APPLYING DEVIATORIC LOADING ============ step= 1000 unbalanced force: 0.007284146090760937 sigma2: -0.0 q= 50000.0 axial deformation (%) 0.0 ~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~ ################## Deviatoric phase is finished and saved successfully ################## ============ RECORD AND PLOT DATA ============ /home/ehsan/yade/install/lib/x86_64-linux-gnu/yade-2019-08-08.git-775ae74/py/yade/plot.py:444: MatplotlibDeprecationWarning: The 'verts' kwarg was deprecated in Matplotlib 3.0 and will be removed in 3.2. Use 'marker' instead. scatter=pylab.scatter(scatterPt[0] if not math.isnan(scatterPt[0]) else 0,scatterPt[1] if not math.isnan(scatterPt[1]) else 0,s=scatterSize,color=line.get_color(),**scatterMarkerKw) ************** END ************** Analysis has been taken for 270.053 seconds or 4.500883333333333 minutes [[ ^L clears screen, ^U kills line. F12 controller, F11 3D view (press "h" in 3D view for help), F10 both, F9 generator, F8 plot. ]] In [1]: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Questions: 1- Is there any problem with getting the following message in the output right after "DEFINING ENGINES"? >> The constructor with a shareWidget is deprecated, use the regular contructor >> instead. 2- Why the volume box calculated by triax is near to 0? >> Box Volume: 6.90392657805617e-310 3- Why in the "APPLYING CONFINING PRESSURE" section the porosity is about 0.8 while I defined it as 0.35? >> porosity= 0.8343868539188172 Thank you, Ehsan -- 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
