Hi Yade friends
I do need run below code with two different engine including (Cundall Contact Law, Hertz Mindlin Law). At the stage of one cundall
Contact must be used and then in Stage two Hertz Mindlin Contact law must be used only. Actually, I do not know how to turn off an
engine at the end of stage one and then turn another one with Hertz contact.
###################################
from mpl_toolkits.axes_grid1 import host_subplot
import matplotlib; matplotlib.rc('axes',grid=True)
import matplotlib.pyplot as plt
from yade import pack
import pylab
from numpy import *
import numpy as np
from math import *
utils.readParamsFromTable(seed=1,num_spheres=10000,compFricDegree =8.0)
from yade.params import table
seed=table.seed
num_spheres=table.num_spheres# number of spheres
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
confiningS=-1e5
young1=10e7
psdSizes,psdCumm=[0.00509,0.00704,0.0101,0.0127,0.0166,0.0222,0.0338],[.322,.376,.452,.533,.632,.726,1]
psdSizesArray=np.array(psdSizes)
psdSizesMeter=psdSizesArray*0.001 #Convert the size of particles to meter
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(0.0001,0.0001,0.0001)
sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=seed)
sp.psd(bins=4000,mass=True)
key='_Youngmadules_'+str(young1)
O.materials.append(FrictMat(young=40e7,poisson=.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=40e7,poisson=.3,frictionAngle=0,density=0,label='frictionless'))
walls=aabbWalls((mn,mx),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)
O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
triax=TriaxialStressController(
internalCompaction=False,
goal1=confiningS,
goal2=confiningS,
goal3=confiningS,
max_vel=10,
label="triax"
)
newton=NewtonIntegrator(damping=0)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_FrictMat_FrictMat_MindlinPhys(betan=0.04,betas=0.04)],
[Law2_ScGeom_FrictPhys_CundallStrack(label = 'Cundall_Law'),Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=True,includeMoment=False,label
= 'Mindlin_Law'),
]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
VTKRecorder(iterPeriod=1000,fileName="./export",recorders=['spheres','intr','color'],dead=1,label='VTK'),
triax,
TwoPhaseFlowEngine(label="unsat"),
newton
]
always_use_Cundall_Law=True
always_use_Mindlin_Law=False
while 1:
O.run(1000,True)
unb=unbalancedForce()
print ('unb=', unb,'meanstress:',abs(triax.goal1-triax.meanStress)/abs(triax.goal1),'e=', triax.porosity/(1-triax.porosity), )
if unb<0.5 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.5 and (triax.porosity)/(1-triax.porosity)<0.979999999999999999 :
break
print ("Particle.Distribution.is.stable")
##################################
##################################
finalFricDegree = 9.0
triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))
while 1:
O.run(1000,True)
unb=unbalancedForce()
if unb<0.5 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.5:
break
ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2]
si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2]
from yade import plot
O.engines=O.engines+[PyRunner(iterPeriod=20,command='history()',label='recorder')]
## a function saving variables
def history():
plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2,
s11=-triax.stress(0)[0]-si0,
s22=-triax.stress(2)[1]-si1,
s33=-triax.stress(4)[2]-si2,
ee=(((triax.porosity)/(1-triax.porosity))*1e-2),
syy=((-triax.stress(3)[1])*1e-3),
i=O.iter
)
#plot.addData( i=O.iter,ee=(triax.porosity)/(1-triax.porosity),s22=-triax.stress(3)[1]*1e-3)
#plot.addData(e22=-triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5)))
##make nice animations:
#O.engines=O.engines+[PyRunner(iterPeriod=200,command='flow.saveVtk()')]
from yade import plot
plot.plots={'syy':('ee')}
plot.plot()
###################################
#Stage one ( Two phase flow )
###################################
print("Get saturation")
triax.stressMask=7
triax.goal1=triax.goal3=0 #confiningS
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=0
void=[];sy=[];Suction=[];Void2=[];valumstrain=[]
meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]
unsat.isPhaseTrapped=True
unsat.initialization()
unsat.surfaceTension = 0.0728
#file=open('Results'+key+'.txt',"w")
for pg in arange(1.0e-5,2,0.7):
unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]
unsat.invasion()
unsat.computeCapillaryForce()
for b in O.bodies:
O.forces.setPermF(b.id, unsat.fluidForce(b.id))
while 1:
O.run(1000,True)
unb=unbalancedForce()
if unb<0.5:
break
print ('PC:',((-unsat.bndCondValue[2])*1e-3), 'Sw:',(unsat.getSaturation(False)),'e:',(triax.porosity)/(1-triax.porosity), 'e33=', (triax.strain[2]), 'e22=', (-triax.strain[1]), 's22=', (-triax.stress(3)[1]*1e-3))
print ('----------LOADING--------------')
###############################################################
###Stage Two ( Loading )
###############################################################
O.engines[2].lawDispatcher.functors[1].always_use_Cundall_Law=False
O.engines[2].lawDispatcher.functors[1].always_use_Mindlin_Law=True
unsat.bndCondIsPressure=[0,0,0,1,0,0]
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.internalCompaction=False
triax.wall_bottom_activated=False
loadingMatrix=[-5e3,-2e5,-8.3e5,-4e5]
for i in arange(0,len(loadingMatrix),1):
triax.goal2=loadingMatrix[i]
O.run(1000,True)
void.append((triax.porosity)/(1-triax.porosity))
sy.append((-triax.stress(3)[1])*1e-3)
valumstrain.append((-triax.strain[1]))
print ('e=',((triax.porosity)/(1-triax.porosity)), 's22=',(-triax.stress(3)[1]*1e-3), 'Load:',(loadingMatrix[i]), 'Sw:',(unsat.getSaturation(False)), 'e22=',(-triax.strain[1]), 'PC=',((-unsat.bndCondValue[2])*1e-3))
plt.figure()
plt.subplot(221)
plt.grid(True)
plt.plot(sy,void,linestyle='-',marker="^", markersize=8,color='r',linewidth=2)
plt.legend(('Exprimental data','simulation data'),loc='upper right', shadow=True)
plt.show()
