#simulation that provides stress, porosity and permeability within all the points of the sample for each flux from yade import pack,export import yade.timing import math from pylab import rand import os.path import numpy as np import matplotlib.pyplot as plt ##### PREAMBLE FOR BATCH EXEC ##### # load parameters from file if run in batch # default values are used if not run from batch readParamsFromTable(normalizedGap=0.8,param2=55.5) #param2 is just an example # make normalizedGap, param2 accessible directly as variables later from yade.params.table import * #For using create a file named e.g. 'param.table' with parameters in columns, typical content: #normalizedGap param2 #0.5 .1 #0.333 .2 #0.2 .3 #then run: yadedaily-batch -j=2 --job-threads=1 param.table injectionPeriodic.py #'-j' defines the total number of cores, '--job-threads' defines the number of core per task (1 task= 1 line of the batch file) ##### END PREAMBLE ##### O.periodic=True #Input values width=0.8 aperture=0.05*width gap = width*normalizedGap height=1 depth=0.4 numSpheres=2000 #Originally 5000 in the script from E.P. Montella compFricDegree=10 rate=0.1 damp=0.0 mu=0.01 key='EduTest'+str(normalizedGap) key2='EduTest2' porosity=0.8 now=-0.001 newSample=False#turn this true if you want to generate new sample by pluviation each time you run the script #Deduced values filename="init"+key+str(numSpheres) volume = width*height*depth meanRad=pow(volume*(1-porosity)/(math.pi*(4/3.)*numSpheres),1/3.) if (meanRad*6>depth): print("INCOMPATIBLE SIZES. INCREASE DEPTH OR INCREASE NUM_SPHERES") if not os.path.isfile(filename+".yade") or newSample:#we create new sample if it does not exist... O.cell.hSize=Matrix3(width, 0, 0, 0 ,3.*height, 0, 0, 0, depth) O.materials.append(FrictMat(young=400000.0,poisson=0.5,frictionAngle=compFricDegree/180.0*math.pi, density=1600,label='spheres')) O.materials.append(FrictMat(young=400000.0,poisson=0.5,frictionAngle=radians(15), density=1000,label='walls')) lowBox = box( center=(width/2.0,height,width/2.0), extents=(width*1000.0,0,width*1000.0) ,fixed=True,wire=False) O.bodies.append(lowBox) topBox = box( center=(width/2.0,2*height+4*meanRad,width/2.0), extents=(width*1000.0,0,width*1000.0) ,fixed=True,wire=False) O.bodies.append(topBox) sp=pack.SpherePack() sp.makeCloud((0,height+2*meanRad,0),(width,2*height+2*meanRad,depth),-1,.002,numSpheres,periodic=True,porosity=porosity,seed=2) O.bodies.append([sphere(s[0],s[1],color=(0.6+0.15*rand(),0.5+0.15*rand(),0.15+0.15*rand()),material='spheres') for s in sp]) else:#... else we re-use the previous one O.load(filename+".yade") Gl1_Sphere.stripes=True O.usesTimeStepper=True flow=PeriodicFlowEngine( dead=1, meshUpdateInterval=40, defTolerance=-1, debug=0, permeabilityFactor=1.0, useSolver=3, duplicateThreshold=depth*6, wallIds=[-1,-1,0,1,-1,-1], bndCondIsPressure=[0,0,0,1,0,0], viscosity=0.01, label="flow" ) ################################################################## ##analytical expression and difference between analytical and DEM solution ##analytical solution anfloat=width #sample width bnfloat=height*0.32 #sample height ###analytical solution from math import * #import numpy #import numpy as np #from math import log #from math import sqrt a=anfloat b=bnfloat delta = 0.04 x = np.arange(0, a+delta, delta) y = np.arange(0, b+delta, delta) X, Y = np.meshgrid(x, y) c=b*2 ####################### effective stress field xmin=0 xmax=width ymin=height ymax=height+height nx=math.trunc(a/delta) ny=math.trunc(b/delta) l=5*meanRad #from yade import utils #from yade import * #from math import * def weight(r,l): return exp(-(abs(r/l)**2)) def averageS(xmin,xmax,ymin,ymax,nx,ny,l): list=[] list2=[] s=utils.bodyStressTensors() TW=PeriodicFlowEngine(duplicateThreshold=depth,wallIds=[-1,-1,0,1,-1,-1]) for i in range(0,nx+1): x=xmin+i*(xmax-xmin)/float(nx) for j in range(0,ny+1): y=ymin+j*(ymax-ymin)/float(ny) meanStress=0; weights=0; volumeTotal=0;volumeSolid=0 for b in O.bodies: if type(b.shape)==Sphere: r= ((x-b.state.pos[0])**2+(y-b.state.pos[1])**2)**0.5 w=weight(r,l) weights+=w volumeSphere=4.*pi/3.*b.shape.radius**3 volumeCell=TW.volume(b.id) #note: next expression is actually w*volumeCell*[volumeSolid/volumeCell], i.e. the equivalent continuum stress meanStress+=w*volumeSphere*s[b.id][1,1] volumeTotal+=w*volumeCell volumeSolid+=w*volumeSphere b.spherePoro=1.-volumeSphere/volumeCell localPorosity=1-volumeSolid/volumeTotal list.append([x,y,meanStress/volumeTotal,localPorosity]) list2.append([x,y,localPorosity]) return list,list2 def plot_porosity(): data = averageS(xmin,xmax,ymin,ymax,nx,ny,l)[1] y = [data[i][1] for i in range(ny+1)] x = [data[i*ny][0] for i in range(nx+1)] poro = [data[i][2] for i in range((nx+1)*(ny+1))] grid_poro = np.zeros((ny+1,nx+1)) for i in range(nx+1): for j in range(ny+1): grid_poro[j,i] = poro[i*(ny+1)+(ny-j)] #print(data) #print(grid_poro) fig = plt.figure(figsize=(20, 10), facecolor='w', edgecolor='w') #plt.imshow(grid_poro, extent=[np.amin(x), np.amax(x), np.amin(y), np.amax(y)], vmin=np.amin(abs(grid_poro)), vmax=np.amax(abs(grid_poro))) plt.imshow(grid_poro, extent=[np.amin(x), np.amax(x), np.amin(y), np.amax(y)], vmin=0.35, vmax=0.65) """interpolation='nearest', aspect='auto', cmap='seismic'""" plt.xlabel('x [m]') plt.ylabel('y [m]') plt.colorbar() plt.show() def spherePoro(xmin,xmax,ymin,ymax,nx,ny,l): s=utils.bodyStressTensors() TW=PeriodicFlowEngine(duplicateThreshold=depth,wallIds=[-1,-1,0,1,-1,-1]) for b in O.bodies: if type(b.shape)==Sphere: volumeSolid=4.*pi/3.*b.shape.radius**3 volumeTotal=TW.volume(b.id) b.spherePoro=1.-volumeSolid/volumeTotal def vitpos(): velo=[] posi=[] hs=0.9*height for i in range(len(O.bodies)): if (O.bodies[i].state.pos[1])>hs: velo.append(abs(O.bodies[i].state.vel[1])) posi.append(O.bodies[i].state.pos[0]) return velo,posi ################################################################## newton=NewtonIntegrator(damping=damp,gravity=(0,-9.81,0)) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],allowBiggerThanPeriod=1), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1000,timestepSafetyCoefficient=0.3, defaultDt=utils.PWaveTimeStep(),label='timestepper'), flow, newton, PyRunner(command="flow.saveVtk(key)",iterPeriod=10,dead=1,label="vtkE"), VTKRecorder(recorders=["spheres","velocity","stress"],iterPeriod=10,dead=1,fileName=key+'_',label="vtkR") ] ## get static equilibrium and save the result for later use if not os.path.isfile(filename+".yade") or newSample: O.run(1000,1) while unbalancedForce()>0.01: O.run(100,1) O.save(filename+".yade") ##################################Permeability############################## ## determining the value of permeability considering Darcy's law and hydraulic head. Hydrostatic pressure is neglected #flow function when null pressure is considered on the top and 200Pa are considered at the bottom of the sample flowPermUp=PeriodicFlowEngine( dead=0, meshUpdateInterval=40, defTolerance=-1, debug=0, permeabilityFactor=1.0, useSolver=3, duplicateThreshold=depth, wallIds=[-1,-1,0,1,-1,-1], bndCondIsPressure=[0,0,1,1,0,0], bndCondValue=[0,0,10,0,0,0], boundaryUseMaxMin=[0,0,0,0,0,0], viscosity=0.01 ) ###insert them as extra engines for one step, just to get permeability ##O.engines=O.engines+[flowPermUp, ##flowPermDown] ##O.run(1,1) ##O.engines=O.engines[:-2] #remove last two ##Qin = flowPermUp.getBoundaryFlux(0) ##Qout = flowPermUp.getBoundaryFlux(1) ##per_up= abs(flowPermUp.getBoundaryFlux(0))*mu/((width*depth))/(10/(height*0.32-2*meanRad)) ##print("Qin_up=",Qin," Qout_up=",Qout," permeability_Up=",per_up) ##Qin = flowPermDown.getBoundaryFlux(0) ##Qout = flowPermDown.getBoundaryFlux(1) ##per_down= abs(flowPermDown.getBoundaryFlux(1))*mu/((width*depth))/(10/(height*0.32-2*meanRad)) ##print("Qin_down=",Qin," Qout_down=",Qout," permeability_down=",per_down) ##per=(per_down+per_up)/2 ##print("Average permeability=",per) ##vtkE.dead=vtkR.dead =0 ########### find particle at the top of the sample ############## #simpleHo=0 #for s in O.bodies: #if isinstance(s.shape,Sphere): #pos=s.state.pos #if pos[1]>simpleHo: #simpleHo=pos[1] ########### define what to plot ############## #from yade import plot #import pylab #def myAddPlotData(): #index = flow.getCell(0.5*width,height,0.5*depth) #if index>0: ########### find particle at the top of the sample ############## #simpleH=0 #for s in O.bodies: #if isinstance(s.shape,Sphere): #pos=s.state.pos #if pos[1]>simpleH: #simpleH=pos[1] ########### function to compute hf ############## #cavityh=height #for s in O.bodies: #v=s.state.vel #magvel=pow((v[0]*v[0]+v[1]*v[1]+v[2]*v[2]),0.5) #if magvel>0.14: #pos=s.state.pos #if pos[1]>cavityh: #cavityh=pos[1] ########################## #plot.addData(t=O.time,i=O.iter, #p=flow.getPorePressure((0.5*width,height,0.5*depth)), #q=-now, #Ho=simpleHo-1, #hf=cavityh-1, #H=simpleH-1) #H=simpleH-1 #hf=cavityh-1 ##remove the extra flow engines (last two positions), add a recorder #O.engines=O.engines[:-2]+[PyRunner(command=('myAddPlotData()'),iterPeriod=50,label="recorder")] #plot.plots={'t':('p',None,('q','b')),'t ':('hf','H',('Ho','r--')),'q':('hf','H',('Ho','r--')),'q ':'p'} #flow.dead=0 #O.saveTmp() ##plot.plot() #commented out in batch #def imposeFlux(valF, gap):#precondition #found=0 #listF1=[] #listF2=[] #for k in range(flow.nCells()): #if 0 in flow.getVertices(k) and flow.getCellCenter(k)[0]>((width-aperture-gap)/2.) and flow.getCellCenter(k)[0]<((width+aperture-gap)/2.): #listF1.append(k) #if 0 in flow.getVertices(k) and flow.getCellCenter(k)[0]>((width-aperture+gap)/2.) and flow.getCellCenter(k)[0]<((width+aperture+gap)/2.): #listF2.append(k) #flow.clearImposedFlux() #if len(listF1)==0: flow.imposeFlux((0.5*width,height,0.5*depth),valF) #elif len(listF2)==0: flow.imposeFlux((0.5*width,height,0.5*depth),valF) #else: #for k in listF1: flow.imposeFlux(flow.getCellBarycenter(k),valF/float(len(listF1))) #for k in listF2: flow.imposeFlux(flow.getCellBarycenter(k),valF/float(len(listF2))) ##print "flux imposed in ",len(listF)," cells" ###########################impose the costant flux######################## #flux=-0.01 #this one is large and it will fluidize violently, see below for smoother hydro-loading #k=1 #now=flux*k #flow.blockHook="imposeFlux(now,gap)" ##VTE=export.VTKExporter("./test") ##spherePoro(xmin,xmax,ymin,ymax,nx,ny,l) ##VTE.exportSpheres(what=[("poro","b.spherePoro")]) ##########################increasing the flux############################### ##flux=-0.0001 ##resPi=[] ##resFi=[] ##resHi=[] ##reshfi=[] ##k=0 ##while k<40: ##now=flux-0.0005*k ##print ('iter=',k) ##flow.blockHook="imposeFlux(now)" ##O.run(2000,True) ##if (k==0): ##stresslist=averageS(xmin,xmax,ymin,ymax,nx,ny,l) ##for i in range(len(stresslist)): ##vector=stresslist[i] ##if (vector[0]==width/2) and (vector[1]==height): ##porosity1 = vector[3] ##pressure=flow.getPorePressure((0.5*width,height,0.5*depth)) ##print ('pressure=',pressure) ##alpha=per*((1-porosity1)**2)/((porosity1)**3) ##print ('alpha=',alpha) ##print ('porosity=',porosity1) ##porosity1=0 ##qflux=2*now*mu/(2*pi*per)/(depth) ##N=100. ##sources=[-aperture/2.+aperture*u/(N-1.)-0.001 for u in range(int(N))] ##B=0 ##for sx in sources: ##A=0 ##for i in range(-L,L+1): ##for j in range(-H,H+2): ##A = A+(log(((sqrt(((0.5*anfloat-(a+0.001)/2-i*a-sx)**2+(0-j*c)**2))))))*(-1)**(abs(j)) ##B=A+B ##pressure_y_max=(qflux/N)*B ##stresslist=averageS(xmin,xmax,ymin,ymax,nx,ny,l) ##file = open('test_stress_0.1width_N5000_increasing.txt', 'a') ##txt = 'x' + '\t' +'y' +'\t' + 'stress' +'\t'+'porosity'+'\t'+'Perm'+'\t'+'flux'+'\n' ##file.write(txt) ##for i in range(len(stresslist)): ##vector=stresslist[i] ##porosity1 = vector[3] ##per3=alpha*(porosity1**3)/((1-porosity1)**2) ##txt = str(vector[0]) + '\t'+ str(vector[1]) + '\t' + str(vector[2]) +'\t' + str(vector[3]) +'\t' + str(per3) +'\t' + str(-now) +'\n' ##file.write(txt) ### Close your file ##file.close() ##k=k+1 ##########################decreasing the flux############################## ##flux=-0.0125 ##resPd=[] ##resFd=[] ##resHd=[] ##reshfd=[] ##k=0 ##while k<26: ##flow.clearImposedFlux() ##now=flux+0.0005*k ##flow.imposeFlux((0.5*width,height,0.5*depth),now) ##O.run(6000 ,True) ##porosity1=0 ##stresslist=averageS(xmin,xmax,ymin,ymax,nx,ny,l) ##file = open('test_stress_0.1width_N5000_decreasing.txt', 'a') ##txt = 'x' + '\t' +'y' +'\t' + 'stress' +'\t'+'porosity'+'\t'+'Perm'+'\t'+'flux'+'\n' ##file.write(txt) ##for i in range(len(stresslist)): ##vector=stresslist[i] ##porosity1 = vector[3] ##per3=alpha*(porosity1**3)/((1-porosity1)**2) ##txt = str(vector[0]) + '\t'+ str(vector[1]) + '\t' + str(vector[2]) +'\t' + str(vector[3]) +'\t' + str(per3) +'\t' + str(-now) +'\n' ##file.write(txt) ### Close your file ##file.close() ##k=k+1 ###plot.saveGnuplot("increasing_decreasing_4") ##################################################################