mindlin-Deresiewicz contact law

Hi,

The contact law you are attempting to use is in fact not available. I started coding it during my PhD but I never completed it since I decided that the simplified HM law was good enough for my purposes. Please see below the link with the list of contact laws currently in use with Yade and which have been tested and accurately verified by means of elements contact tests.
https://www.yade-dem.org/wiki/ConstitutiveLaws

Please feel free to add your contact law should you need the complete path described by the HM law (and not the simplified version which is now in the code).

Cheers,
Chiara

On 27 Dec 2013, at 18:56, Sina Jafari wrote:

> Question #241108 on Yade changed:
> https://answers.launchpad.net/yade/+question/241108
>
> Status: Answered => Open
>
> Sina Jafari is still having a problem:
> for example in the following script I have modified the triaxial example script a little bit to study the effect of contact law on the output, the point is when I use mindlin-Deresiewicz contact law, no shear force is recorded in the output which is in contradiction with the results obtained by, say, cundall law. Why is this happening? do I need to define a parameter or somethin'? Thanks for your answers. here is the script:
> # -*- coding: utf-8 -*-
> #*************************************************************************
> # Copyright (C) 2010 by Bruno Chareyre *
> # bruno.chareyre_at_grenoble-inp.fr *
> # *
> # This program is free software; it is licensed under the terms of the *
> # GNU General Public License v2 or later. See file LICENSE for details. *
> #*************************************************************************/
>
> ## This script details the simulation of a triaxial test on sphere packings using Yade
> ## See the associated pdf file for detailed exercises
> ## the algorithms presented here have been used in published papers, namely:
> ## * Chareyre et al. 2002 (http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)
> ## * Chareyre and Villard 2005 (https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
> ## * Scholtès et al. 2009 (http://dx.doi.org/10.1016/j.ijengsci.2008.07.002)
> ## * Tong et al.2012 (http://dx.doi.org/10.2516/ogst/2012032)
> ##
> ## Most of the ideas were actually developped during my PhD.
> ## If you want to know more on micro-macro relations evaluated by triaxial simulations
> ## AND if you can read some french, it is here: http://tel.archives-ouvertes.fr/docs/00/48/68/07/PDF/Thesis.pdf
>
> from yade import pack,plot
> import matplotlib; matplotlib.rc('axes',grid=True)
> import pylab
> ############################################
> ### DEFINING VARIABLES AND MATERIALS ###
> ############################################
> key='_Kenney_'
> num_spheres=48710
> psdSizes,psdCumm=[0.262*0.89185,0.53*0.89185,0.97*0.89185,1.76*0.89185,2.49*0.89185,3.4*0.89185,4.87*0.89185,6.4*0.89185],[0.1,3.4,9.1,19.1,29.4,48.7,81.5,100]
> #targetPorosity = 0.387 #the porosity we want for the packing
> compFricDegree = 26.5 # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
> finalFricDegree = 26.5 # contact friction during the deviatoric loading
> rate=0.0001 # loading rate (strain rate)
> damp=0.2 # damping coefficient!!!!!!!!!!
> stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
> young=540e6 # contact stiffness
> mn,mx=Vector3(0,0,0),Vector3(53.52,53.52,53.52) # corners of the initial packing
>
>
> ## create materials for spheres and plates
> O.materials.append(FrictMat(young=young,poisson=0.35,frictionAngle=radians(compFricDegree),density=2000,label='spheres'))
> O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
>
> # create walls around the packing
> walls=aabbWalls([mn,mx],material='walls',oversizeFactor=1)
> wallIds=O.bodies.append(walls)
>
> ## use a SpherePack object to generate a random loose particles packing
> sp=pack.SpherePack()
> sp.particleSD2(radii=psdSizes,passing=psdCumm,numSph=7500,cloudPorosity=0.55)
> O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
> #or alternatively (higher level function doing exactly the same):
> #sp.toSimulation(material='spheres')
>
> ############################
> ### DEFINING ENGINES ###
> ############################
>
> triax=TriaxialStressController(
> ## ThreeDTriaxialEngine will be used to control stress and strain. It controls particles size and plates positions.
> ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
> maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)!!!!!!
> finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)!!!!!!!!!!
> thickness = 0,
> ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
> ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
> ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
> ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
> ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
> 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(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_MindlinPhys(frictAngle=26.5)],
> [Law2_ScGeom_MindlinPhys_MindlinDeresiewitz()]
> ),
> ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
> triax,
> TriaxialStateRecorder(iterPeriod=100,file='WallStresses.dat'),
> newton
> ]
>
> ########################################
> #### APPLYING CONFINING PRESSURE ###
> ########################################
>
> #the value of (isotropic) confining stress defines the target stress to be applied in all three directions
> 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.01:
> break
>
> O.save('confinedState'+'.yade.gz')
> print "### Isotropic state saved ###"
> print 'ACN=',utils.avgNumInteractions(),'Porosity=',utils.voxelPorosityTriaxial(triax),'Calculation Time(Sec)=',O.realtime
>
> ####################################################
> #### REACHING A SPECIFIED POROSITY PRECISELY ###
> ####################################################
>
> ### We will reach a prescribed value of porosity with the REFD algorithm
> ### (see http://dx.doi.org/10.2516/ogst/2012032 and
> ### http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)
>
> #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(500,1)
>
> #O.save('compactedState'+key+'.yade.gz')
> #print "### Compacted state saved ###"
> #print 'ACN=',utils.avgNumInteractions(),'Porosity=',utils.voxelPorosityTriaxial(triax),'Calculation Time(Sec)=',O.realtime
>
> ##############################
> ### 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 = 0
> #now goal2 is the target strain rate
> triax.goal2=-rate
> # we define three lateral stresses during the test, here the same 10kPa as for the initial confinement.
> triax.goal1=-rate
> triax.goal3=-rate
>
> #we can change damping here. What is the effect in your opinion?
> #newton.damping=0.1
>
> #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 avgcoord():
> stress=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
> plot.addData(ACN=utils.avgNumInteractions(),
> P=stress,e=utils.voxelPorosityTriaxial(triax)/(1-utils.voxelPorosityTriaxial(triax)))
> plot.saveDataTxt('results C-MAO')
> def addplotdirect():
> utils.plotDirections(noShow=True).savefig('Interaction histogram.jpeg')
>
> f = open('%s.txt'%"interactions of each body",'w')
> def numint():
> i=O.iter
> f.write('%s\n'%' ')
> f.write('%s\n'%' ')
> f.write('%s\n'%' ')
> f.write('%s\n'%int(i))
> f.write('%s\n'%float(triax.meanStress))
> f.write('%s\t'%str('bodyid'))
> f.write('%s\t'%str(" "))
> f.write('%s\t'%str('Radius'))
> f.write('%s\t'%str(" "))
> f.write('%s\n'%str('NumInteractions'))
> for m in O.bodies:
> intrs=m.intrs()
> nintrs=len(intrs)
> bodyid=m.id
> if isinstance(m.shape,Sphere):
> radii=m.shape.radius
> else:
> radii="NaN"
> f.write('%s\t'%int(bodyid))
> f.write('%s\t'%str(" "))
> f.write('%s\t'%float(radii))
> f.write('%s\t'%str(" "))
> f.write('%s\n'%int(nintrs))
>
> g = open('%s.xls'%"Contact forces",'w')
> def cntctforce():
> stress=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
> if (stress>0.995*200000 and stress<1.005*200000) or (stress>0.995*400000 and stress<1.005*400000):
> g.write('%s\n'%int(O.interactions.countReal()))
> g.write('%s\n'%float(triax.meanStress))
> for j in O.interactions:
> if not j.isReal: continue
> fn = j.phys.normalForce.norm()
> fs = j.phys.shearForce.norm()
> g.write('%s\t'%float(fn))
> g.write('%s\n'%float(fs))
>
>
>
>
> # include a periodic engine calling that function in the simulation loop
> O.engines=O.engines[0:5]+[PyRunner(iterPeriod=100,command='cntctforce()')]+[PyRunner(iterPeriod=50000,command='addplotdirect()')]+[PyRunner(iterPeriod=50000,command='numint()')]+[PyRunner(iterPeriod=5000,command='avgcoord()')]+O.engines[5:10]
> #plot.plots={'P':('ACN',None,'e')}
> #plot.plot()
>
> O.run(100,True)
> O.run(1000000)
> #### PLAY THE SIMULATION HERE WITH "PLAY" BUTTON OR WITH THE COMMAND O.run(N) #####
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp