New question #663499 on Yade: https://answers.launchpad.net/yade/+question/663499
Hi all, Back to my previous question #663085 - I have updated the script as per below - it is working without any error however when I change the percentage of materials from 10% to 90%( through this line :- if random.random() < 0.9:) the stress- strain results does not change at all. Even I am not sure if "DEVIATORIC LOADING" phase is impaling, because the stress value does not go more than 100kpa in all different ratios of the materials combinations. Would you please advise me where is the issue? Thanks Seti #### #!/usr/bin/python # -*- coding: utf-8 -*- from __future__ import division from yade import plot,pack,timing import time, sys, os, copy from yade import ymport, utils, pack, export #from pylab import * import math import numpy as np ############################################ ### DEFINING VARIABLES AND MATERIALS ### ############################################ # default parameters or from table readParamsFromTable(noTableOk=True, # unknownOk=True, targetPorosity=0.45, #the porosity we want for the packing num_spheres=3000,# number of spheres young=236e9, poisson=.2, density=2600, frictionAngle=radians(30), normalCohesion=3.e6, shearCohesion=1.e6, etaRoll=0.1, compFricDegree=30, # contact friction during the confining phase finalFricDegree=30, # contact friction during the deviatoric loading key='_triax_base_', # put you simulation's name here rate=-0.02, # loading rate (strain rate) damp=0.7, # damping coefficient stabilityThreshold=0.01, # we test unbalancedForce against this value in different loops (see below) ) from yade.params import table from yade.params.table import * mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing ## create materials for spheres and plates mat1=CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,normalCohesion=normalCohesion,shearCohesion=shearCohesion,momentRotationLaw=True,etaRoll=etaRoll,label='cement') mat_1=O.materials.append(mat1) mat2=CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,normalCohesion=3.e3,shearCohesion=1.e3,momentRotationLaw=True,etaRoll=etaRoll,label='glass_sphere') mat_2=O.materials.append(mat2) O.materials.append(FrictMat(young=2*young,poisson=.25,frictionAngle=0,density=0,label='frictionlessWalls')) ## create walls around the packing walls=aabbWalls([mn,mx],thickness=0,material='frictionlessWalls') wallIds=O.bodies.append(walls) ## use a SpherePack object to generate a random loose particles packing sp=pack.SpherePack() sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same O.bodies.append([sphere(center,rad,material='cement') for center,rad in sp]) cement=[] glass_sphere=[] for b in O.bodies: if not isinstance(b.shape,Sphere): # change material only on spheres continue if random.random() < 0.9: b.mat = mat1 b.shape.color = (1,0,0) b.state.mass*=mat1.density/mat1.density cement.append(b.id) else: b.mat = mat2 b.shape.color = (0,1,1) b.state.mass*=mat2.density/mat1.density glass_sphere.append(b.id) ############################ ### DEFINING ENGINES ### ############################ triax=TriaxialStressController( maxMultiplier=1.+2.4e5/young, # spheres growing factor (fast growth) finalMaxMultiplier=1.+2.4e4/young, # spheres growing factor (slow growth) thickness = 0, stressMask = 7, internalCompaction=True, # If true the confining pressure is generated by growing particles ) newton=NewtonIntegrator(damping=damp) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( #box-sphere interactions will be the simple normal-shear law, we use ScGeom for them [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()], #Boxes will be frictional (FrictMat), so the sphere-box physics is FrictMat vs. CohFrictMat, the Ip type will be found via the inheritance tree (CohFrictMat is a FrictMat) and will result in FrictPhys interaction physics #and will result in a FrictPhys [Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp"), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(normalCohesion=MatchMaker(matches=((mat_1,mat_1,3.e6),(mat_1,mat_2,3.e4),(mat_2,mat_2,3.e3))),shearCohesion=MatchMaker(matches=((mat_1,mat_1,1.e6),(mat_1,mat_2,1.e4),(mat_2,mat_2,1.e3))))], #Finally, two different contact laws for sphere-box and sphere-sphere [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment( useIncrementalForm=True, #useIncrementalForm is turned on as we want plasticity on the contact moments always_use_moment_law=False, #if we want "rolling" friction even if the contact is not cohesive (or cohesion is broken), we will have to turn this true somewhere label='cohesiveLaw')], ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), TriaxialStateRecorder(iterPeriod=1000,file='90%,matchmaker,j'+key), triax, newton ] triax.goal1=triax.goal2=triax.goal3=-100000 while 1: O.run(1000, True) ##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium unb=unbalancedForce() print 'unbalanced force:',unb,' mean stress: ',triax.meanStress if unb<stabilityThreshold and abs(-100000-triax.meanStress)/100000<0.001: break print 'timestep at confinedState:',O.iter import sys #this is only for the flush() below while triax.porosity>targetPorosity: ## we decrease friction value and apply it to all the bodies and contacts compFricDegree = 0.95*compFricDegree setContactFriction(radians(compFricDegree)) #print "\r Friction: ",compFricDegree," porosity:",triax.porosity, sys.stdout.flush() ## while we run steps, triax will tend to grow particles as the packing ## keeps shrinking as a consequence of decreasing friction. Consequently ## porosity will decrease O.run(1000,1) print 'porosity:',triax.porosity print 'timestep at compactedState:',O.iter print "### Compacted state saved ###" ############################## ### DEVIATORIC LOADING ### ############################## ##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant triax.internalCompaction=False ## Change contact friction (remember that decreasing it would generate instantaneous instabilities) setContactFriction(radians(finalFricDegree)) ##set stress control on x and z, we will impose strain rate on y triax.stressMask = 5 ##now goal2 is the target strain rate triax.goal2=rate ## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement. triax.goal1=-100000 triax.goal3=-100000 newton.damping=0.1 O.engines[2].lawDispatcher.functors[1].always_use_moment_law = True #We assign cohesion to all contacts at the next iteration O.engines[2].physDispatcher.functors[1].setCohesionNow = True O.run(1000,True) ##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button. O.saveTmp() ##################################################### ### Example of how to record and plot data ### ##################################################### #from yade import plot ### a function saving variables 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) if 1: ## include a periodic engine calling that function in the simulation loop O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7] ##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder')) else: ## With the line above, we are recording some variables twice. We could in fact replace the previous ## TriaxialRecorder ## by our periodic engine. Uncomment the following line: O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder') O.run(100,True) ### declare what is to plot. "None" is for separating y and y2 axis #plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')} ### the traditional triaxial curves would be more like this: #plot.plots={'e22':('s11','s22','s33',None,'ev')} plot.plots={'e22':('s11','s22')} ## display on the screen (doesn't work on VMware image it seems) plot.plot() ##### PLAY THE SIMULATION HERE WITH "PLAY" BUTTON OR WITH THE COMMAND O.run(N) ##### ## In that case we can still save the data to a text file at the the end of the simulation, with: plot.saveDataTxt('90%,matchmaker '+key) ##or even generate a script for gnuplot. Open another terminal and type "gnuplot plotScriptKEY.gnuplot: plot.saveGnuplot('90%,matchmaker'+key) rr=yade.qt.Renderer() rr.shape=True rr.intrPhys=True -- 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