New question #663066 on Yade: https://answers.launchpad.net/yade/+question/663066
Dear all, I have done some small changes on cohesive triaxial script - changing the normal and cohesion values for two type of materials - but the results does not make sense - unbalace force nan and mean stress zero can you please advise me what is my mistake thanks #!/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.e6,shearCohesion=1.e6,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.5: 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")], #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), 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) 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(1000,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('400,alphaKr=0.6,normalCohesion=1e9,shearCohesion=.3e9,young=236e9,j'+key) ##or even generate a script for gnuplot. Open another terminal and type "gnuplot plotScriptKEY.gnuplot: plot.saveGnuplot('400,alphaKr=0.6,normalCohesion=1e9,shearCohesion=.3e9,young=236e9,j'+key) rr=yade.qt.Renderer() rr.shape=True rr.intrPhys=True Thanks -- 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