New question #688653 on Yade: https://answers.launchpad.net/yade/+question/688653
Hello all, I am using Ubuntu 18.04, and Yade 2019-08-08.git-775ae74 I developed a Triaxial test based on [1] in the following form. This code is working with rigid walls/boundaries which move in the "deviatoric" part of the test and this leads to the deformation of the specimen. In this case, I will not be able to capture the shear band and localization. To do so, I need to define the walls/boundaries as the deformable/flexible boundaries. In fact, I should not let the boundaries to move, but I should let them be deformed based on their flexibility. My question is that what do I need to do to make these rigid walls to flexible boundaries? [1] https://gitlab.com/yade-dev/trunk/blob/master/examples/triax-tutorial/script-session1.py My code: print ('************** START **************') import numpy as np import time import datetime, os start_time=time.time() from datetime import datetime import math from yade import qt, export, utils from yade import pack ################################## ######### DEFINING VARIABLES ######### nRead=readParamsFromTable( num_spheres=20000, compFricDegree = 29, key='_triax_', unknownOk=True ) from yade.params import table num_spheres=table.num_spheres key=table.key targetPorosity = 0.4 compFricDegree = table.compFricDegree finalFricDegree = 29 IP=100 # iteration period to record data and stuff micro_record_iterPeriod=IP ORN=3000 # O.Run Number of iterations micro_record_enable_normal_branch=True micro_record_float_type = np.float32 damp=0.2 thick=0 stabilityThreshold=0.01 PCPC=0.0001 # Precision of Confining Pressure Convergence r_min=0.1*1e-3 # m d_min=2*r_min # m r_max=0.3*1e-3 # m d_max=2*r_max # m r_avr=(r_min+r_max)/2 # m d_avr=2*r_avr # m r_fuz=(r_max/r_avr)-1 # m Kn=10e8*(d_avr) ### FIXME Kt=10e8*(d_avr) ### FIXME young=Kn/r_avr # 2 (E r1 E r2 / E r1 + E r2) >>> E = Kn/r_avr poisson=Kn/Kt # Kt/Kn Ls=0.02 # m length of specimen ### FIXME L_REV=7*(d_avr) # m if Ls < L_REV: sys.exit("*** ERROR! The specimen's dimension is too samll! ***") elif Ls==L_REV: print ("*** This is the minimum specimen's dimension you can take! ***") else: print ("*** The specimen's dimension is good enough! ***") mn,mx=Vector3(0,0,0),Vector3(Ls,Ls,Ls) Vt=-1*1e-3 # m/s # negative sign describes the compression direction strainRate=Vt/Ls # %/sec target_strain=0.25 ### FIXME % print ("The target strain has been set to:", target_strain) sigmaIso=-5e5 # Pa ### FIXME particleDensity=2000 #kg/m3 ################################################## ################# DEFINING MATERIALS ################# O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=particleDensity,label='spheres')) O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls')) ################################################ ################# DEFINING PACKING ################# walls=aabbWalls([mn,mx],thickness=thick,material='walls') for w in walls:w.shape.radius=0 wallIds=O.bodies.append(walls) sp=pack.SpherePack() clumps=False sp.makeCloud(mn,mx,r_avr,r_fuz,num_spheres,False, 0.95,seed=1) O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp]) from yade import export os.mkdir('3axresults') os.chdir('/home/ehsan/Desktop/3axresults') export.text('InitialPackingData') #################################################### ################# DEFINING TRIAXIAL TEST ################# triax=TriaxialStressController( maxMultiplier=1.+2e4/young, finalMaxMultiplier=1.+2e3/young, thickness = thick, stressMask = 7, internalCompaction=True, ) ################################################## ################# 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, t = O.time, # virtual (yade) time --- time of simulation fab = utils.fabricTensor()[0]) ################################################ ################# DEFINING ENGINES ################# 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()] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), triax, PyRunner(iterPeriod=IP,command='history()',label='macro_recorder'), TriaxialStateRecorder(iterPeriod=IP,file='WallStresses'+table.key), newton ] Gl1_Sphere.stripes=True if nRead==0: yade.qt.Controller(), yade.qt.View() ########################################################## ################# APPLYING CONFINING PRESSURE ################# triax.goal1 = sigmaIso triax.goal2 = sigmaIso triax.goal3 = sigmaIso while 1: O.run(ORN,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 ConfStressRatio=abs(sigmaIso-triax.meanStress)/abs(sigmaIso) print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~') print ('mean stress engine:',triax.meanStress,' mean stress (Calculated):',meanS, ' ConfSratio:',ConfStressRatio,' step:', O.iter/ORN, ' Time:',O.time, ' TimeStep',O.dt) print ('porosity:',triax.porosity, ' void ratio:',triax.porosity/(1-triax.porosity)) if unb<stabilityThreshold and ConfStressRatio<PCPC: break export.text('FinalPhase01PackingData') e22Check=-triax.strain[1] ###%%%*** print ('Axial Strain',e22Check) print ('Mean stress engine: ',triax.meanStress) print ('Mean stress (Calculated): ',meanS) ######################################################## ################# REACHING TARGET POROSITY ################## print ('============ REACHING TARGET POROSITY ============') import sys while triax.porosity>targetPorosity: compFricDegree = 0.95*compFricDegree setContactFriction(radians(compFricDegree)) print ('\r Friction: ',compFricDegree,' porosity:',triax.porosity,'step= ',O.iter/ORN,' Time:',O.time, ' TimeStep',O.dt) sys.stdout.flush() O.run(ORN,True) ################################################## ################# DEVIATORIC LOADING ################# triax.internalCompaction=False setContactFriction(radians(finalFricDegree)) triax.wall_top_activated=True triax.wall_bottom_activated=False triax.wall_left_activated=True triax.wall_right_activated=True triax.wall_back_activated=True triax.wall_front_activated=True triax.stressControl_1=True triax.stressControl_2=False triax.stressControl_3=True triax.strainRate2=strainRate triax.stressMask = 5 triax.goal1 = sigmaIso triax.goal2 = strainRate triax.goal3 = sigmaIso newton.damping=0.1 ### FIXME os.mkdir('Phase02PackingData') os.chdir('/home/ehsan/Desktop/3axresults/Phase02PackingData') while 1: O.run(ORN,True) export.text(str(O.iter/ORN)) record_micro_data(str(O.iter/ORN)) unb=unbalancedForce() axialS=triax.stress(triax.wall_top_id)[1] eps2=triax.strain[1] eps2cal=(triax.height-Ls)/Ls print ('~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~') print ('sigma2: ',axialS, ' q: ', axialS-sigmaIso,' step= ', O.iter/ORN,' Time:',O.time, ' TimeStep',O.dt) print ('Strain2 Calculated: ',eps2cal, ' Strain2 engine: ',eps2,'axial deformation (%)', (eps2-e22Check)*100) print('A:', abs(eps2-e22Check), ' B:',target_strain) if abs(eps2-e22Check)>=target_strain: #triax.sigma2 break ###################################### ################# END ################# print ('************** END **************') O.realtime print ('Analysis has been taken for',O.realtime, 'seconds or', O.realtime/60, 'minutes') print('Real time of run:',(time.time() - start_time), 'seconds or',(time.time() - start_time)/60, 'minutes') yade.timing.stats() ####################################################### ################# RECORD Macro DATA and Plot ################# print ('============ RECORD AND PLOT DATA ============') O.run(ORN,True) plot.plots={'e22':('ev',),'j':('s22',),'e22 ':('s22',),' j':('s11','s33',)} plot.labels={'s11':'$\sigma_{11}$' , 's22':'$\sigma_{22}$' , 's33':'$\sigma_{33}$' , 'e22':'$\epsilon_{22}$' , 'ev':'$\epsilon_{V}$' , 'i':'Time Step'} plot.plot() os.chdir('/home/ehsan/Desktop/3axresults') plot.saveDataTxt('Macro_results') O.pause() %%%%%%%%%%%%%%%%%%%%% Thank you for your helps. -- 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
