yade-users team mailing list archive

[ehsan benabbas <question688165@xxxxxxxxxxxxxxxxxxxxx>]

Question #688165 on Yade changed:
https://answers.launchpad.net/yade/+question/688165

ehsan benabbas gave more information on the question:
This is my code:

######################################################################################################
######### TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z IS THE FRONT AXIS #########
######################################################################################################

print ('************** START **************')
import numpy as np
import datetime
import time
import math
from yade import qt, export, utils
from datetime import datetime

######################################
######### DEFINING VARIABLES #########

print ('============ DEFINING VARIABLES ============')

nRead=readParamsFromTable(
 num_spheres=20000,
 compFricDegree = 30,
 key='_triax_base_',
 unknownOk=True
)

from yade.params import table

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.4
compFricDegree = table.compFricDegree
finalFricDegree = 30
rate=-0.02
damp=0.2
thick=0
stabilityThreshold=0.01
young=200000000
poisson=1
mn,mx=Vector3(0,0,0),Vector3(0.02,0.02,0.02)
sigmaIso=-500000
particleDensity=2000
######################################################
################# DEFINING FUNCTIONS #################

print ('============ 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,
    n = triax.porosity)

######################################################
################# DEFINING MATERIALS #################

print ('============ 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 #################

print ('============ DEFINING PACKING ============')
walls=aabbWalls([mn,mx],thickness=thick,material='walls')
wallIds=O.bodies.append(walls)

from yade import pack
sp=pack.SpherePack()
clumps=False
volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
sp.makeCloud(mn,mx,0.0002,0.5,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

from yade import export
export.text('PackingData')

##########################################################
################# DEFINING TRIAXIAL TEST #################

print ('============ DEFINING TRIAXIAL TEST ============')
triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = thick,
 stressMask = 7,
 internalCompaction=True,
)

####################################################
################# DEFINING ENGINES #################

print ('============ DEFINING ENGINES ============')
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=1,command='history()',label='recorder'),
 TriaxialStateRecorder(iterPeriod=1,file='WallStresses'+table.key),
 NewtonIntegrator(damping=damp,label="newton")
]
Gl1_Sphere.stripes=True
if nRead==0: yade.qt.Controller(), yade.qt.View()
print ('Number of elements: ', len(O.bodies))
print ('Box Volume: ', triax.boxVolume)
print ('Box Volume calculated: ', volume)

###############################################################
################# APPLYING CONFINING PRESSURE #################

print ('============ APPLYING CONFINING PRESSURE ============')
triax.stressmask=7
triax.goal1 = sigmaIso
triax.goal2 = sigmaIso
triax.goal3 = sigmaIso
while 1:
    O.run(1000, 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
    print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~')
    print ('unbalanced force:',unb,' mean stress engine: ',triax.meanStress,' mean stress (Calculated): ',meanS)
    print ('porosity=',triax.porosity)
    print ('void ratio=',triax.porosity/(1-triax.porosity))
    if unb<stabilityThreshold and abs(sigmaIso-triax.meanStress)/sigmaIso<0.001:
        break
O.save('confinedPhase'+key+'.xml')
e22Check=triax.strain[1]
print ('Axial Strain',e22Check)
print ('Mean stress engine: ',triax.meanStress)
print ('Mean stress (Calculated): ',meanS)
print ('################## Isotropic phase is finished and saved successfully ##################')

#############################################################
################# REACHING TARGET POROSITY ##################

print ('============ REACHING TARGET POROSITY ============')
import sys #this is only for the flush() below

while triax.porosity>targetPorosity:
	compFricDegree = 0.95*compFricDegree
	setContactFriction(radians(compFricDegree))
	print ("\r Friction: ",compFricDegree," porosity:",triax.porosity)
	sys.stdout.flush()
	O.run(1000,1)

print ('Final Number of elements: ', len(O.bodies))
print ('Final Box Volume: ', triax.boxVolume)
print ('Final Box Volume calculated: ', volume)
print ('Final porosity=',triax.porosity)
print ('Final void ratio=',triax.porosity/(1-triax.porosity))
O.save('compactedState'+key+'.yade.gz')
print ('################## Target porosity is reached and compacted state saved successfully ##################')

######################################################
################# DEVIATORIC LOADING #################

print ('============ APPLYING DEVIATORIC LOADING ============')
triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))
triax.stressControl_2 = 0
triax.stressMask = 5
triax.goal1 = sigmaIso
triax.goal2 = rate
triax.goal3 = sigmaIso
newton.damping = 0.1
unb = unbalancedForce()
axialS = triax.stress(triax.wall_top_id)[1]
print ('step=', O.iter, 'unbalanced force:',unb,' sigma2: ',axialS, 'q=', axialS-sigmaIso)
print ('Axial Deformation (%)', (triax.strain[1]-e22Check)*100)
print ('~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~')
O.save ('final.xml')
print ('################## Deviatoric phase is finished and saved successfully ##################')

########################################################
################# RECORD AND PLOT DATA #################

print ('============ RECORD AND PLOT DATA ============')
if 1:
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=1,command='history()',label='recorder')]+O.engines[5:7]
else:
  O.engines[4]=PyRunner(iterPeriod=1,command='history()',label='recorder')
O.run(100,True)
plot.plots={'e22':('ev')}
plot.labels={'s11':'$\sigma_{11}$' , 's22':'$\sigma_{22}$' , 's33':'$\sigma_{33}$' , 'e22':'$\epsilon_{22}$' , 'ev':'$\epsilon_{V}$'}
plot.plot()
plot.saveDataTxt('results'+key)
plot.saveGnuplot('plotScript'+key)

#####################################################
################# RECORD Micro DATA #################

data = []
for i in O.interactions:
   fn = i.phys.normalForce
   fs = i.phys.shearForce
   cp = i.geom.contactPoint
   normal = i.geom.normal
   b1,b2 = [O.bodies[id] for id in (i.id1,i.id2)]
   p1,p2 = [b.state.pos for b in (b1,b2)]
   branch = p2 - p1
   cp,normal,branch,fn,fs = [tuple(v) for v in (cp,normal,branch,fn,fs)] # Vector3 -> tuple
   d = dict(cp=cp,normal=normal,branch=branch,fn=fn,fs=fs)
   data.append(d)
# new data contains the information, you can save it e.g. as JSON
import json
with open("interactions.json","w") as f:
   json.dump(data,f,indent=3)

# =============================================================================
# with open(fileName,"w") as ff:
#    ff.write("# cpx cpy cpz fnx fny fnz ...\n")
#    for d in data: # or directly for i in O.interactions
#       cp = d["cp"]
#       fn = d["fn"]
#       ...
#       ff.write("{} {} {} {} {} {} ... "\n".format(cp[0],cp[1],cp[2],fn[0],fn[1],fn[2],...)
# =============================================================================
#######################################
################# END #################

print ('************** END **************')
O.realtime
print ('Analysis has been taken for',O.realtime, 'seconds or', O.realtime/60, 'minutes')


Cheers,
Ehsan

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.