New question #691645 on Yade:
https://answers.launchpad.net/yade/+question/691645
I’m working on a Fluid-DEM coupling simulation script and trying to export two
things : (1) the total stress acting at the boundary and (2) the effective
stress with changing pore water pressure.
What the script will do :
State I : Execute an iso-compression using PeriTriaxController engine
State II : The PeriTriaxController engine is switched off, the
PeriodicFlowEngine is then applied with fluid pressure increasing by stepping
value Delta_P = 100 kPa after every 1000 iterations. Once the fluid pressure
reaches 500 kPa , the simulation stops and a stress plot will be given.
The script runs fine without any technical error ; however, when running with
YADE 2018.02b (installed from Ubuntu repo), the output stress is calculated as
it’s supposed to be : the boundary stress is constant while the effective
stress decreases according to the imposed fluid pressure; but when running with
Yadedaily (version Yade 20200701-4022~5cf8771~bionic1 at spoken time) , I
receive a different behaviour: the boundary stress increases gradually.
Can you please tell me what changes have been made in YADE between the two
versions that produces such different results on the same script? What
modification should be done in latest Yade version to obtain the result given
in YADE 2018.2b?
More info:
1) I’m running Ubuntu 18.04.04
2) useSolver was set to zero because Yade 2018.2b gives me error about
not-compiled CHOLMOD when I used 3, but I think I should ask in another
question about this issue.
The MWE script is below:
from __future__ import print_function
from yade import pack,plot,export,utils
import math
sp=pack.SpherePack()
O.periodic=True
# dimensions of sample (fixed by particle size such as L/D~X)
DIAMETER =1.e-1
RADIUS =0.5*DIAMETER
length =10*(DIAMETER)
height =length
width =length
thickness =length/100.
# microproperties
DENS =2600
E =1.e9
P =0.25
compFRIC =1.
FRIC =30.
TENS =0.
COH =0.
# boundary conditions
PI =1.e5 # for sample preparation: isotropic compaction up to PI pressure
PC =PI # confining pressure for assembly preparation
SN =5.e6 # normal stress
# simulation control
DAMPSHEAR =0.
ITER =1e6
OUT ='test'
#### create sample and loading boxes
O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)
O.materials.append(CohFrictMat(isCohesive=True,density=DENS,young=E,poisson=P,frictionAngle=radians(0.),normalCohesion=1e100,shearCohesion=1e100,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=DENS,young=E,poisson=P,frictionAngle=radians(compFRIC),normalCohesion=TENS,shearCohesion=COH,label='sphereMat'))
upBox =
utils.box(center=(length/2.0,2*height+thickness/2.0,width/2.0),orientation=Quaternion(1,0,0,0),extents=(length,thickness/2.,width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox =
utils.box(center=(length/2.0,height-thickness/2.0,width/2.0),orientation=Quaternion(1,0,0,0),extents=(length,thickness/2.,width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
O.bodies.append([upBox,lowBox])
sp.makeCloud((0,height+1.5*RADIUS,0),(length,2*height-1.5*RADIUS,width),rMean=RADIUS,rRelFuzz=0.2,periodic=True)
O.bodies.append([utils.sphere(s[0],s[1],color=(0,0.5,1),material='sphereMat')
for s in sp])
effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol
#### engines
flow=PeriodicFlowEngine(
isActivated=0 # no flow calculation if this is False
,useSolver=0 # value 3 not work with YADE 2018.2b package version
,defTolerance=+1 # with this value, triangulation is done according to
meshUpdateInterval
,meshUpdateInterval=1000
,duplicateThreshold=0.5
,boundaryUseMaxMin=[0,0,1,1,0,0]
,wallIds=[-1,-1,1,0,-1,-1]
,wallThickness=thickness
,bndCondIsPressure=[0,0,0,0,0,0]
,bndCondValue=[0,0,0,0,0,0]
,permeabilityFactor=1
,viscosity=1e-3
,fluidBulkModulus=2.2e9
,label='flowEng'
)
O.engines=[
ForceResetter()
,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
,InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
)
,flow
,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep(),label='timeStepper')
,NewtonIntegrator(damping=0.3,label='newton')
,PyRunner(command='dataRecorder()',iterPeriod=1,label='recData',dead=True)
]
inputP=0
def dataRecorder():
global inputP
h=vol=vol_s=nb_s=0.
h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
contactStress=getStress(vol)
for o in O.bodies:
if isinstance(o.shape,Sphere):
nb_s += 1
vol_s += 4.*pi/3.*(o.shape.radius)**3
n = 1-vol_s/vol
nbFrictCont=0.
for i in O.interactions:
if i.isReal and i.phys.cohesionBroken:
nbFrictCont+=1
plot.addData(
iter=O.iter
,stress_upWall0=abs(O.forces.f(0)[0]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall1=abs(O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])),stress_upWall2=abs(O.forces.f(0)[2]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
,contactStress00=(contactStress[0,0]),contactStress01=(contactStress[0,1]),contactStress02=(contactStress[0,2]),contactStress10=(contactStress[1,0]),contactStress11=abs(contactStress[1,1]),contactStress12=(contactStress[1,2]),contactStress20=(contactStress[2,0]),contactStress21=(contactStress[2,1]),contactStress22=(contactStress[2,2])
,x=O.bodies[0].state.pos[0]
,pf=inputP
,height=h
,volume=vol
,porosity=n
,k=2.0*nbFrictCont/nb_s
,Ek=kineticEnergy()
,unbF=unbalancedForce()
)
phase=0
def triaxDone():
global phase
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])
if phase==0:
h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
contactStress=getStress(vol)
vol_s=Rmean=Rmax=nbSph=0
Rmin=1e6
for o in O.bodies:
if isinstance(o.shape,Sphere):
nbSph+=1
Rmean+=o.shape.radius
if o.shape.radius>Rmax: Rmax=o.shape.radius
if o.shape.radius<Rmin: Rmin=o.shape.radius
vol_s += 4.*pi/3.*(o.shape.radius)**3
Rmean=Rmean/nbSph
n = 1-vol_s/vol
print('DONE! iter =',O.iter,'| sample generated: nb spheres =',nbSph,',
Rmean =',Rmean,', Rratio =',Rmax/Rmin,', porosity =',n)
print('Changing contact properties now')
utils.setContactFriction(radians(FRIC))
print('APPLYING CONFINING PRESSURE: sx,sy and sz will go to PC =',PC)
triax.goal=(-PC/volRatio,-PC/volRatio,-PC/volRatio) # - - -
phase+=1
elif phase==1:
print('DONE! iter =',O.iter,'| isotropic confinement done: stresses
=',volRatio*triax.stress)
triax.dead=True
O.pause()
#### Initialization
print('SAMPLE PREPARATION!')
O.run(1000000,1)
O.step()
#### Applying normal stress
print('NORMAL LOADING! iter =',O.iter)
stage=0
stiff=fnPlaten=currentSN=0.
def servo():
global stage,stiff,fnPlaten,currentSN
if stage==0:
currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
unbF=unbalancedForce()
boundaryVel=copysign(min(0.1,abs(0.5*(currentSN-SN))),currentSN-SN)
O.bodies[0].state.vel[1]=boundaryVel
if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
stage+=1
fnPlaten=O.forces.f(0)[1]
print('Normal stress =',currentSN,' | unbF =',unbF)
for i in O.interactions.withBody(O.bodies[0].id):
stiff+=i.phys.kn
print('Normal stiffness =',stiff)
print('DONE! iter =',O.iter)
O.pause()
if stage==1:
fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
boundaryVel=copysign(min(0.1,abs(0.35*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired)
# why 0.35?
O.bodies[0].state.vel[1]=boundaryVel
O.engines =
O.engines[:4]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[4:]
O.run(1000000,1)
recData.dead=False
O.run(1,1)
#### injecting fluid
print('FLUID! iter =',O.iter)
flowEng.isActivated=1
flowEng.boundaryUseMaxMin=[1,1,0,0,1,1]
flowEng.bndCondIsPressure=[0,0,1,0,0,0]
newton.damping=0.
O.run(1,1)
# automatic
iter0=O.iter
iterMax=5e3
deltaP=1e5
while O.iter<(iter0+iterMax):
O.run(1000,True)
inputP+=deltaP
flowEng.bndCondValue=[0.,0.,+inputP,0.,0.,0.]
flowEng.updateBCs()
print('updateBCs! inputP=',inputP)
O.run(1,1)
plot.plots={'iter':('pf','stress_upWall1','contactStress11')}
plot.plot()
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
