New question #693568 on Yade: https://answers.launchpad.net/yade/+question/693568
Hello everyone, I'm trying to simulate a periodic undrained triaxial test by mantaining constant volume instead of explicitly simulating water (as a sidenote, I never could get periodicFlowEngine to work correctly for this). My simulation is divided in three phases, for the first one a cloud of spheres is isotropically compacted to a target porosity, then for the second phase the packing is isotropically compacted again to a target stress, and in the third one a velGrad is imposed in the cell so there is compression a certain rate in one direction and there is extension at half the compression rate in the other two directions. The idea of this final phase is to impose a constant volume condition as to simulate an undrained behaviour. The problem is that even though at the end of the isoCompaction phase the stress is equal to a target stress (100000 for expample) and volume in the shear phase is indeed constant, after a number of iterations stresses drop to zero and start to oscillate around that value. This is the result I get from my MWE that is at the end of this post: Welcome to Yade 20201018-4279~61d08ab~bionic1 Using python version: 3.6.9 (default, Oct 8 2020, 12:12:24) [GCC 8.4.0] TCP python prompt on localhost:9000, auth cookie `acusye' XMLRPC info provider on http://localhost:21000 qt5ct: using qt5ct plugin Running script /home/geotesis/DEM/Yade/PeriodicModel_noWater_MWE.py Start: 2020-10-21 16:09:53.985841 Cloud created, number of particles = 480 Compacting to target porosity: 2020-10-21 16:09:54.003525 [[ ^L clears screen, ^U kills line. F12 controller, F11 3D view (press "h" in 3D view for help), F10 both, F9 generator, F8 plot. ]] Target porosity passed, porosity = 0.49986744040494263 Starting isocompaction simulation: 2020-10-21 16:10:16.938606 Current iteration = 1534900 , s11 = -100157.84795970531 , s22 = -100145.52886621605 , s33 = -100199.74978550954 , e = 0.6462475530919375 Starting shear simulation: 2020-10-21 16:12:28.666616 Current iteration = 1534901 , s11 = -100154.30036085597 , s22 = -100148.62230105326 , s33 = -100200.53394271889 , e = 0.6462475530919375 Current iteration = 1834901 , s11 = -55393.437466646115 , s22 = -65305.28482971604 , s33 = -51735.14392560943 , e = 0.6462749983880726 Current iteration = 2134901 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , e = 0.6462749983315953 Current iteration = 2434901 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , e = 0.6462749981327709 Current iteration = 2734901 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , e = 0.6462749980919655 Current iteration = 3034901 , s11 = -0.11991207951660474 , s22 = -0.22419400597402317 , s33 = -0.026624167711505737 , e = 0.6462749980486365 Current iteration = 3334901 , s11 = -0.047911199264998226 , s22 = -0.00023733397986893054 , s33 = -0.013137608101300035 , e = 0.6462749979832385 Current iteration = 3634901 , s11 = -0.018439277204848392 , s22 = -0.026673245415610117 , s33 = -0.049101414208267816 , e = 0.6462749979320861 Current iteration = 3934901 , s11 = -0.00045092480566355026 , s22 = -0.0041975100217730744 , s33 = -0.07083510893219645 , e = 0.6462749978859719 Current iteration = 4234901 , s11 = -10.470723417543306 , s22 = -5.878484222842999 , s33 = -2.042145920074655 , e = 0.6462749977701003 Current iteration = 4534901 , s11 = -0.014448912716733427 , s22 = -0.07253827703395103 , s33 = -0.013568831616481059 , e = 0.6462749976951889 Current iteration = 4834901 , s11 = -0.004931408116028016 , s22 = -0.03483611806022092 , s33 = -0.0062370493144523085 , e = 0.6462749976127734 Current iteration = 5134901 , s11 = 0.00012374680500377676 , s22 = -0.0009835150202924818 , s33 = -0.0049291412638700184 , e = 0.6462749975535781 Current iteration = 5434901 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , e = 0.6462749974492507 Finished simulation phase: 2020-10-21 16:18:16.761407 As you can see, after 600000 iterations stresses in all directions drop from 100000 to zero. When I visually check what's happening in the simulation, It's quite obvious that particles start to separate in this phase and start to float around until the target strains are reached. I know that there's something I'm not doing or that I'm doing wrong, but I can't wrap my head around the fact that volume stays constant but particles have space to float around at the same time. This also happens if the isoCompaction goal is higher, though it takes more iterations. As a summary, my CU TX periodic simulation has zero stresses in shear phase, even though it started from isotropic confinement. Expected behaviour would be that stresses change from initial confinement (s11 and s33 equal, and s22 different), wihout all of them reaching zero, while volume stays constant. # -*- coding: utf-8 -*- ###################### ### INITIALIZE ### ###################### from yade import pack from yade import plot ## Time of simulation import time import datetime print ("Start: ",datetime.datetime.now()) start = time.time() ## These are default parameters for batch execution in full code confinement=100 e0=1.0 rate=0.01 batch = 0 ####################################### ### BASIC SIMULATION PARAMETERS ### ####################################### ### Set Periodic cell O.periodic=True x1 = 0.006 y1 = x1 z1 = x1 mn,mx = Vector3(0,0,0),Vector3(x1,y1,z1) # corners of the initial packing O.cell.setBox((x1,y1,z1)) ### PSD in [mm] # Silica Sand #7 psdSizes = [0.075,0.106,0.250,0.425,0.850] psdCumm = [3,7,88,100,100] ## From [mm] to [m] and scale in homotetic way psdScale = 3 # 1 is around 10000 particles, 2 around 1500, 3 around 500 psdSizes = [x*(10**-3)*psdScale for x in psdSizes] psdCumm = [x*(10**-2) for x in psdCumm] ## Material of spheres young=70e9 compFricDegree = 3 # initial contact friction during the confining phase finalFricDegree = 30 # contact friction during the deviatoric loading O.materials.append(FrictMat(young=young, poisson=0.17, frictionAngle=radians(compFricDegree), density=2600, label="spheres")) ################### ### ENGINES ### ################### sigmaIso = -1e8 triax = PeriTriaxController( # specify target values and whether they are strains or stresses goal = (sigmaIso,sigmaIso,sigmaIso), stressMask = 7, # type of servo-control dynCell = True, mass = 0.5, maxStrainRate = (1,1,1), # wait until the unbalanced force goes below this value maxUnbalanced = 0.01, relStressTol = 0.002, doneHook="changePhase()" ) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_MindlinPhys()], [Law2_ScGeom_MindlinPhys_Mindlin()] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,label="GSTS"), triax, NewtonIntegrator(damping = 0.2), PyRunner(command = "checkPorosity()", iterPeriod = 1000, label = "checkPorosityPy"), PyRunner(iterPeriod=10000,command="checkGoals()", label="goalChecker",dead=1) ] ### Functions ### phase = 1 [e11, e22, e33] = [0, 0, 0] [e110, e220, e330] = [0, 0, 0] def changePhase(): # Phase 1: Isotropic compaction to reach target porosity # Phase 2: Isotropic compaction to reach confinement # Phase 3: Triax shear global phase global saveL global confinement global rate global key global sigmaIso if phase == 1: # End targetPorosity phase print ("Target porosity passed, porosity = ", utils.porosity()) checkPorosityPy.dead = 1 # Begin isoCompaction phase prepareIsoCompaction(confinement) elif phase == 2: ## Start Shear phase phase = 3 # Add engine to show that stresses go to zero O.engines = O.engines+[PyRunner(iterPeriod=300000, command="showValues()", label="values")] prepareShear(rate,confinement) elif phase == 3: stopSimulation() else: print("There's something wrong with value of phase, ", phase) O.pause() def checkPorosity(): if utils.porosity()<targetPorosity: changePhase() def prepareIsoCompaction(confinement): global phase # Begin isoCompaction phase phase = 2 setContactFriction(radians(finalFricDegree)) O.cell.setBox(O.cell.size) triax.maxStrainRate = (0.1,0.1,0.1) sigmaIso = -1e3*confinement # kPa triax.goal[0]=triax.goal[1]=triax.goal[2]=sigmaIso print("Starting isocompaction simulation: ", datetime.datetime.now()) def prepareShear(rate,confinement): global phase global e110,e220,e330 # Save current strain [e110, e220, e330] = triax.strain # So things move in the correct direction triax.stressMask = 7 triax.goal[1] = -1e3*confinement triax.goal[0] = triax.goal[2] = -triax.goal[1]/2 # Impose strainRate (du/dx,du/dy,du/dz,dv/dx,dv/dy,dv/dz,dw/dx,dw/dy,dw/dz) O.cell.velGrad=Matrix3(rate/2,0,0,0,-rate,0,0,0,rate/2) triax.maxStrainRate = [rate/2,rate,rate/2] goalChecker.dead = 0 showValues() print("Starting shear simulation: ", datetime.datetime.now()) def showValues(): # This function is just for the MWE print("Current iteration = ", O.iter, ", s11 = ", triax.stress[0], ", s22 = ", triax.stress[1], ", s33 = ", triax.stress[2], ", e = ", utils.porosity()/(1-utils.porosity())) def checkGoals(): global phase global batch if phase == 3: # These values are just for the MWE chk1 = triax.strain[0] - e110 > 0.025 chk2 = triax.strain[1] - e220 < -0.05 chk3 = triax.strain[2] - e330 > 0.025 if chk1 & chk2 & chk3: stopSimulation() else: print("There's something wrong with value of phase, ", phase) O.pause() def stopSimulation(): print("Finished simulation phase: ", datetime.datetime.now()) O.pause() ################################# ### START SIMULATION ### ################################# ### Create initial cloud with desired PSD (distributed by mass) sp = SpherePack() sp.makeCloud(mn,mx, psdSizes=psdSizes, psdCumm=psdCumm, periodic=True, distributeMass=True, seed=1) print ("Cloud created, number of particles = ", len(sp)) ## Send spheres to simulation sp.toSimulation(material="spheres") ### Reach target porosity ### print ("Compacting to target porosity: ", datetime.datetime.now()) ## Porosity of cloud and targetPorosity targetVoidRatio = e0 targetPorosity = (targetVoidRatio)/(1+targetVoidRatio) ### Start simulation ### O.dt = utils.PWaveTimeStep() O.run() -- 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
