Question #556258 on Yade changed:
https://answers.launchpad.net/yade/+question/556258
Tanvir Hossain posted a new comment:
Hi Jan, here is the script.
#New version, amelioration wrt 20150715 :
#- IMPORTANT Modify the length of the vector evaluated from ndimz to ndimz.
Works with Yade6 and is not adapted to the other cases
#- Add a limiter to avoid negative tau transmitted to the fluid (unphysical and
source of numerical instabilities
#- Add a new configuration for the turbulentFluctuation function. Correspond to
a new formulation of the C++ HydroForceEngine. Possibility to use a new
turbulent fluctuation model where each particle fluctuation lifetime depends on
its wall-normal position.
#- Add an option for fluidProfile imposed case.
#- Register the averaged velocity of the part1 and 2 in twoSize mode (Yade6)
#If activated, adapt the script to the associated computer path
Cluster = 0
philippeComputer = 0
julienComputer = 0
lamiaComputer = 0
raphaelComputer = 0
#Import libraries
from yade import pack, plot
import math
import random as rand
import numpy as np
np.set_printoptions(threshold=numpy.nan)
import yade.timing
import pylab as pyl
from matplotlib import pyplot as pyp
import matplotlib.gridspec as gridspec
#Register the path for later use (tests + save)
scriptPath = os.path.abspath(os.path.dirname(sys.argv[-1]))
##
## Main parameters of the simulation
##
diameterPart = 6e-3 #Diameter of the particles, in meter
densPart = 2500 #density of the particles, in kg/m3
waterDepth = 7. #Water depth in diameter
ndimz = 525 #Number of cells in the height
from nsmp1d_yade3D525Phi061 import *
lengthCell = 15 #Streamwise length of the periodic cell
widthCell = 15 #Spanwise length of the periodic cell
Nlayer = 10. #nb of layer of particle, in diameter
fluidHeight = (Nlayer+waterDepth)*diameterPart #Height of the flow from
the bottom of the sample
quasi2D = 0
twoSize = 0
#If twoSize option activated, overwrite the previous input #tanvir
if twoSize:
Nlayer1 = 10 #nb of layer of particle type 1
Nlayer2 = 1 #nb of layer of particle type 2
diameterPart1 = 6e-3 #Diameter of the particles type 1, in meter
diameterPart2 = 4e-3 #Diameter of the particles type 2, in meter
fluidHeight = Nlayer1*diameterPart1 + Nlayer2*diameterPart2 +
7*diameterPart1 #free surface position
diameterPart = diameterPart1 #Size 1 taken as reference
#Fluid library importation and precisions
phiPartMax = 0.61 #Value of the maximum volume fraction used for
the thresholding
##
## Option of the simulation
##
#Imposed fluid profile? If yes, require to add a file fluid.py in the folder
with inside the value of vxFluid vector + turbStress vector if you want to add
turbulent fluctuations also.
imposedFluidPro = 0
#Temporary imposed profile: impose a fluid profile for 50 seconds (considered
as equilibrium) and switch to the fluid resolution
initImposedFluidPro = 0
#Activation and choice of the fluid turbulent fluid velocity fluctuation.
Cannot choose both. If none, will not apply fluct.
turbModel1 = 1
turbModel2 = 0
#Stress averaging: if activated, compute and save the stress tensor depth
profile at each execution of measure()
stressSave = 0
pythonStressCalculation = 0 #use the python function instead of C++
#Option to put lateral walls, to break the biperiodicity and study the
influence of walls.
lateralWalls = 0
#Option for time averaging of the DEM data before fluid resolution
timeAveraging = 0 #1 activated, 0 not activated
nbAv = 5. #number of averaging between the fluid resolution
#Averaging parameter : if activated, evaluate the average using python function
partConcVelocity() instead of the C++ one
pythonAverage = 0
#debug stress mode, record all the quantities necessary to study the momentum
balance in detail.
debugStress = 0
#Video
video = 1 #1, make a movie
if video:
snapshotEngine = 0 #To make the movie later using mencoder
(require the qt.viewer open, obtain what you see)
paraview = 1 #To make the movie later using paraview (do not require
the qt.viewer open, faster)
fps = 1. #frame per second #tanvir
##
## Secondary parameters of the simulation (less changed)
##
#Particles parameters
restitCoef = 0.5 #Restitution coefficient of the particles
partFrictAngle = atan(0.4) #friction angle of the particles, in radian
slope = 0.1 #Angle of the slope in radian #tanvir
# Fluid parameters
expoRZ = 3.1 # Richardson Zaki exponent for the hindrance function
of the drag
dpdx = 0.0 #pressure gradient along streamwise direction
densFluid = 1000. #Density of the fluid
kinematicViscoFluid = 1e-6 #kinematic viscosity of the fluid
# Fluid simulation parameters (RANS), see fluidResolution function
dsig = np.zeros(ndimz)
dz = fluidHeight/(1.0*(ndimz-1))
sig=pyl.linspace(0e0,1e0,ndimz)
for i in range(ndimz-1): # calcul du pas d'espace
dsig[i] = sig[i+1] - sig[i]
#In order to fasten the script, some correspondance are supposed. Should be
respected otherwise it does not work well
dtFluid = 1e-5 #Time step for the fluid resolution
fluidResolPeriod = 1e-2 #1/fluidResolPeriod = frequency of the calculation of
the fluid profile
measurePeriod = (1e-1)*5 # We perform measurement every measurePeriod
seconds #tanvir
savePeriod = 1e1 # We save the simulation (to be able to reload) every
savePeriod seconds.
#Quasi 2D option
if quasi2D ==1:
print '\nQuasi 2D option activated\n'
lateralWalls = 1
widthCell = 6.5/6.
phiPartMax = 0.51
print 'phiPartMax = 0.51\n'
#
## Tests on the script
#
#Width check and Value of the fluid library and size of the system
if widthCell<2 and phiPartMax>0.6:
if widthCell<=1:
print 'Impossible to run simulation with width exactly equal to
the diameter particle size !\n'
exit()
else:
print '\nCareful !! Uncoherence between the 2D character and
the value of phiPartMax !! \n'
#Value of ndimz and dz
if ndimz != int(round(30*fluidHeight/diameterPart)):
print('\n!!Problem dz = d/{0} !!\nUse ndimz = {1} instead to have dz =
d/30 (and change accordingly the fluid library
import)\n'.format(diameterPart/fluidHeight*ndimz,int(round(30*fluidHeight/diameterPart))))
# if quasi2D!=1:
# exit()
#Option compatibility
if measurePeriod<fluidResolPeriod:
print '\nNot possible to have measurePeriod<fluidResolPeriod in the
present configuration, modify the script to make measurement inside measure()
instead of fluidResolution() !\n'
exit()
turbModel = 0
if turbModel1 ==1 or turbModel2==1:
turbModel = 1
if turbModel1 ==1 and turbModel2==1:
print "It is not possible to have the two fluctuations model in
the same time ! \n Modify turbModel1 or turbModel2"
exit()
#Reloading simulation
reloadingSimulation = 0
if os.path.exists(scriptPath +'/paramSave.py'):
reloadingSimulation = 0#1 Out for now, useless, it is not working
##
## Initialization of the parameters
##
# Averaging/Save
vxPart = np.zeros(ndimz) # streamwise particle velocity depth profile
vyPart = np.zeros(ndimz) # spanmwise particle velocity depth profile
vzPart = np.zeros(ndimz) # normal particle velocity depth profile
qsMean = 0 #Mean bead flux
dragTerm = np.zeros(ndimz) # drag term transmitted to the fluid vector
averageDrag = np.zeros(ndimz) # drag field vector
stressProfile =
getStressProfile(1.0,ndimz,dz,0.0,np.zeros(ndimz),np.zeros(ndimz),np.zeros(ndimz))
#Matrix full of zero
positionPartX = np.zeros(len(O.bodies)) #Defined global in order to use
positionPartZ = np.zeros(len(O.bodies)) #Defined global in order to use
velocityPartX = np.zeros(len(O.bodies)) # the saving function
velocityPartZ = np.zeros(len(O.bodies)) # the saving function
if twoSize:
phiPart1 = np.zeros(ndimz) # Concentration vector
vxPart1 = np.zeros(ndimz) # streamwise particle velocity depth profile
vyPart1 = np.zeros(ndimz) # spanwise particle velocity depth profile
vzPart1 = np.zeros(ndimz) # wall-normal particle velocity depth profile
phiPart2 = np.zeros(ndimz) # Concentration vector
vxPart2 = np.zeros(ndimz) # streamwise particle velocity depth profile
vyPart2 = np.zeros(ndimz) # spanwise particle velocity depth profile
vzPart2 = np.zeros(ndimz) # wall-normal particle velocity depth profile
#Initialization of the variables necessary for timeAveraging option
if timeAveraging:
n_dt = 1.
vxPart_dt = np.zeros(ndimz)
phiPart_dt = np.zeros(ndimz)
averageDrag_dt = np.zeros(ndimz)
if twoSize:
phiPart1_dt = np.zeros(ndimz)
averageDrag1_dt = np.zeros(ndimz)
vxPart1_dt = np.zeros(ndimz)
vyPart1_dt = np.zeros(ndimz)
vzPart1_dt = np.zeros(ndimz)
phiPart2_dt = np.zeros(ndimz)
averageDrag2_dt = np.zeros(ndimz)
vxPart2_dt = np.zeros(ndimz)
vyPart2_dt = np.zeros(ndimz)
vzPart2_dt = np.zeros(ndimz)
#Parameter initialized only if it is the initial simulation, otherwise they are
charged
if reloadingSimulation !=1:
#Parameters which should be loaded and not initiallized if we recover
the simulation
phiPart = np.zeros(ndimz) # Concentration vector
turbulentViscosity = np.zeros(ndimz) # Turbulent viscocity
waterDepth = 0.0 #Water depth : height of the "clear" water
vxFluid = np.zeros(ndimz) # Fluid velocity at time step n
turbStress = np.zeros(ndimz) #Turbulent stress vector
fileNumber = 0 #counter to save the files
###########################
## PARAMETERS DECLARATION##
###########################
#Geometrical configuration, define useful quantities
height = 5*fluidHeight #heigth of the periodic cell, bigger than the real
height to take into account jumps
length = lengthCell*diameterPart #length of the stream
width = widthCell*diameterPart #width of the stream
gravityVector = Vector3(9.81*sin(slope),0.0,-9.81*cos(slope))
gravityVectorApplied = Vector3(0.0,0.0,-9.81*cos(slope))
groundPosition = height/4.0
VLayer = dz*length*width #Volume of a layer
partVolume = pi/6.*pow(diameterPart,3)
if twoSize:
partVolume1 = partVolume
partVolume2 = pi/6.*pow(diameterPart2,3)
#Particles contact law parameters
maxPressure =
(densPart-densFluid)*phiPartMax*Nlayer*diameterPart*abs(gravityVector[2])
#Estimated max particle pressure ("hydrostatic")
normalStiffness = maxPressure*diameterPart*1e4 #Evaluate the minimal normal
stiffness to be in the rigid particle limit (cf Roux and Combe 2002)
youngMod = normalStiffness/diameterPart #Young modulus of the particles from
the stiffness wanted. Need that in order to use the restitution coefficient
definition + more practical when polydisperse sample.
poissonRatio = 0.5 #poisson's ratio of the particles. Classical values,
not much influence (see ref in Da Cruz et al 2005)
if twoSize:
youngMod1 = normalStiffness/diameterPart1
youngMod2 = normalStiffness/diameterPart2
O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod1,
poisson=poissonRatio, density=densPart, frictionAngle=partFrictAngle,
label='Mat1'))
O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod2,
poisson=poissonRatio, density=densPart, frictionAngle=partFrictAngle,
label='Mat2'))
#Material
O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod,
poisson=poissonRatio, density=densPart, frictionAngle=partFrictAngle,
label='Mat'))
# if it is the initial simulation create the framework and the simulation loop,
otherwise load it
if reloadingSimulation !=1:
##########################################
## FRAMEWORK CREATION (sphere, wall,...)##
##########################################
#Definition of the semi-periodic cell
O.periodic = True
O.cell.setBox(length,width,height)
#To be compatible with lateral walls
leftLimitY = 0.0
rightLimitY = width
centerLimitY = width/2.0
#If lateralWalls, add the walls and increase the cell size to avoid
particle touching the walls to appear on the other side
if lateralWalls:
#Warn the user
print '\nlateralWalls option activated: mono-periodic boundary
condition only !\n'
#Increase the cell size
O.cell.setBox(length,width+2*diameterPart,height)
#Modify accordingly the position of the center of the cell and
the wall right and left position
leftLimitY = diameterPart
rightLimitY = width+diameterPart
centerLimitY = diameterPart + width/2.0
#Define the wall and add to the simulation
sidePlaneL = box(center=
(length/2.0,leftLimitY,height/2.0),extents=(2000,0,height*10),fixed=True,wire =
True,color = (1,0,0), material = 'Mat')
sidePlaneR = box(center=
(length/2.0,rightLimitY,height/2.0),extents=(2000,0,height*10.0),fixed=True,wire=True,
material = 'Mat',color = (0,0,1))
O.bodies.append([sidePlaneR,sidePlaneL])
if quasi2D == 1:
# Regular arrangement of spheres sticked at the bottom with
random height
a = range(0,int(length/(diameterPart))) #The length is divided
in diameter of particle
for b in a: #loop creating a set of sphere sticked at the
bottom with a (uniform) random altitude comprised between 0.5 (diameter/12) and
5.5mm (11diameter/12) with steps of 0.5mm. The repartition is made around
groundPosition. Exactly the same set up as the one from the experiment.
n = rand.randrange(0,12,1)/12.0*diameterPart
#Define a number between 0 and 11/12 diameter with steps of 1/12 diameter
(0.5mm in the experiment)
O.bodies.append(sphere((b*diameterPart,centerLimitY,groundPosition -
11*diameterPart/12.0/2.0 + n),diameterPart/2.,color=(0,0,0),fixed =
True,material = 'Mat'))
else:
#Execute the file containing the database which is a 2D list
containing the positions and radius of the particles composing the random fixed
bed. Each position x,y (integer) contain a list of particle with position and
radius at this place
if Cluster==1:
pathFixedBot = '/home/maurin/'
elif philippeComputer==1:
pathFixedBot = '/home/philippe/fixedBot/'
elif tanvir==1:
pathFixedBot = '/home/tanvir/output/'
elif julienComputer==1:
pathFixedBot = '/home/julien/Yade/scriptsRaphael/'
elif lamiaComputer==1:
pathFixedBot = '/home/lamia/periodicStream/'
elif raphaelComputer:
pathFixedBot =
'/home/raphael/these/scripts/twoWayCoupling/'
else:
print '\nNo precised computer used ! Cannot resolve the
fixed bottom file path accordingly ! Precise the computer used at line 12\n'
exit()
if diameterPart==6e-3:
if lengthCell>= 10:
execfile(pathFixedBot+'fixedBotDataLongMONO.py')
else:
execfile(pathFixedBot+'fixedBotDataShort.py')
print '\n Use the short database for
random fixed bottom generation ! \n Be careful, polydisperse bottom !\n'
elif diameterPart==3e-3:
execfile(pathFixedBot+'fixedBotDataMONOd3.py')
elif diameterPart==12e-3:
execfile(pathFixedBot+'fixedBotDataMONOd12.py')
else:
print '\nno fixed bottom of the right size !! quit the
simulation ! \n'
exit()
xMax = len(fixedBottomData) # Maximum size of the data base
along x axis
yMax = len(fixedBottomData[0]) # Maximum size of the data base
along y axis
x = rand.randrange(0,xMax,1) # Generate a uniform random
number between 0 and the size of the data base along x (in diameterPart)
y = rand.randrange(0,yMax,1) # Generate a uniform random
number between 0 and the size of the data base along y (in diameterPart)
#Starting from the random x and y position, loop over the
length and width of the require fixed bed increasing x and y position (in
diameterPart). Create then a sphere starting from 0 for x and y.
for i in range(0,lengthCell):
for j in range(0,int(widthCell)):
for k in
fixedBottomData[(i+x)%xMax][(j+y)%yMax]:
O.bodies.append(sphere(center =
((i+k[0]/diameterPart-int(k[0]/diameterPart))*diameterPart,(j+k[1]/diameterPart-int(k[1]/diameterPart))*diameterPart+leftLimitY,groundPosition+k[2]),radius
= k[3],material = 'Mat',fixed = True,color = (0,0,0)))
# Ground reference and Wall on the side
lowPlane = box(center=
(length/2.0,centerLimitY,groundPosition),extents=(200,200,0),fixed=True,wire=False,color
= (0.,1.,0.),material = 'Mat') #Build a plane to have a reference for the eyes
WaterSurface = box(center=
(length/2.0,centerLimitY,groundPosition+fluidHeight),extents=(2000,width/2.0,0),fixed=True,wire=False,color
= (0,0,1),material = 'Mat',mask = 0) #Build a plane to have a reference for
the eyes of the position of the water surface. No interaction with the sphere
O.bodies.append([lowPlane,WaterSurface]) #add to simulation
#Count the initial number of particles (for later check)
initNumber = len(O.bodies)
#Create a loose cloud of particle inside the cell
partCloud = pack.SpherePack()
if twoSize:
diameterPart = diameterPart1
partVolume = partVolume1
Nlayer = Nlayer1
#Define the number of particle from the density of particle per unit
length
if quasi2D==1:#to keep the correspondance with previous work where I
putter NbPart = Nl*L/d
partNumber = int(Nlayer*lengthCell)
else:
partNumber =
int(Nlayer*phiPartMax*diameterPart*length*width/partVolume) #Volume of beads
#Define the deposition height considering that the packing realised by
make cloud is 0.2
depositionHeight = Nlayer*phiPartMax/0.2*diameterPart #Consider that
the packing realised by make cloud is 0.2
#Create a cloud of partNumber particles
partCloud.makeCloud(minCorner=(0,leftLimitY,groundPosition+diameterPart+1e-4),maxCorner=(length,rightLimitY,groundPosition+depositionHeight),rRelFuzz=0.,
rMean=diameterPart/2.0, num = partNumber)
#Send this packing to simulation with material Mat
partCloud.toSimulation(material='Mat')
partNumber += initNumber
#Same for a second size of particles
if twoSize:
partCloud2 = pack.SpherePack()
if quasi2D==1:#to keep the correspondance with previous work
where I putter NbPart = Nl*L/d
partNumber2 =
int(Nlayer2*lengthCell*diameterPart1/diameterPart2)
else:
partNumber2 =
int(Nlayer2*phiPartMax*diameterPart2*length*width/partVolume2) #Volume of beads
#Define the deposition height considering that the packing
realised by make cloud is 0.2
depositionHeight2 =Nlayer2*diameterPart2*phiPartMax/0.2
partCloud2.makeCloud(minCorner=(0,leftLimitY,groundPosition+depositionHeight),maxCorner=(length,rightLimitY,groundPosition+depositionHeight+depositionHeight2),rRelFuzz=0.,
rMean=diameterPart2/2.0, num = partNumber2)
partCloud2.toSimulation(material='Mat2',color = (0,0,1)) #send
this packing to simulation with material Mat2
#Increment the counters for check and deposition time input
depositionHeight+=depositionHeight2
partNumber+=partNumber2
#Evaluate the deposition time considering the free-fall time of the
highest particle to the ground (overestimation)
depoTime = sqrt(depositionHeight*2/abs(gravityVector[2]))
#Check
if len(O.bodies)<partNumber:
print '\nProblem : too low number of beads in the sample
compared with what was asked !\n'
exit()
# Collect the ids of the spheres which are dynamic and to which we
should add the hydroforce
dv =1e-2 #Scale of the velocity imposed to the particle, no
influence if not extreme case en = 1
idApplyForce = []
for b in O.bodies:
if isinstance(b.shape,Sphere) and b.dynamic:
if widthCell<2: #For the quasi 2D case, need to push
slightly the spheres randomly in order to break the 2D configuration
deltaV = - dv + rand.uniform(0,2*dv)
b.state.vel[1]+=deltaV
idApplyForce+=[b.id]
#########################
#### SIMULATION LOOP#####
#########################
O.engines = [
# Reset the forces
ForceResetter(),
# Detect the potential contacts
InsertionSortCollider([Bo1_Sphere_Aabb(),
Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod
= True),
# Calculate the different interactions
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
,label = 'interactionLoop'),
#Apply an hydrodynamic force to the particles
HydroForceEngine(densFluid = densFluid,viscoDyn =
kinematicViscoFluid*densFluid,zRef = groundPosition,gravity =
gravityVectorApplied,deltaZ = dz,expoRZ = expoRZ,lift = False,nCell =
ndimz,vCell = length*width*dz ,vxFluid = vxFluid,phiPart = phiPart,vxPart =
vxPart,ids = idApplyForce, label = 'addFluidForce', dead = True),
#Time averaging of the DEM data for fluid resolution
PyRunner(command = 'timeAveragingDEMdata()', virtPeriod =
fluidResolPeriod/nbAv, label = 'timeAve', dead = True),
#Calculate the fluid profile using the routine
PyRunner(command = 'fluidResolution(fluidResolPeriod)',
virtPeriod = fluidResolPeriod, label = 'fluidResol', dead = True),
#Calculate the turbulent fluctuations from a DRW random walk
model
PyRunner(command = 'turbulentFluctuation()', virtPeriod = 0.1,
label = 'turbFluct', dead = True),
#Measurement, output files
PyRunner(command = 'measure()', virtPeriod = measurePeriod,
label = 'measurement', dead = True),
# Check if the packing is stabilized, if yes activate the hydro
force on the grains and the slope.
PyRunner(command='gravityDeposition(depoTime)',virtPeriod =
0.01,label = 'packing'),
#GlobalStiffnessTimeStepper, determine the time step
GlobalStiffnessTimeStepper(defaultDt = 1e-4, viscEl =
True,timestepSafetyCoefficient = 0.7, label = 'GSTS'),
# Integrate the equation and calculate the new
position/velocities...
NewtonIntegrator(damping=0.2, gravity=gravityVector,
label='newton')
]
#save the initial configuration to be able to recharge the simulation
starting configuration easily
O.saveTmp()
#run
O.run()
#Add engines to make a video
if video:
if snapshotEngine:
print '\n video snapshot activated : takes snapshot of
the qt.view as you see it\n'
O.pause()
O.engines +=
[SnapshotEngine(virtPeriod=1.0/fps,fileBase=scriptPath+'/figure/',label='snapshooter')]
elif paraview:
print '\n video paraview activated : register vtk file
for later use of paraview (video...)\n'
O.engines +=
[VTKRecorder(virtPeriod=1/fps,recorders=['spheres','colors'],fileName=scriptPath+'/Paraview/p1-')]
else:
print '\nBe careful the video is not activated : need
to specify snapshotEngine or paraview value.\n'
O.pause()
#Add a pyRunner to re-activate the fluctuations at 50s
if initImposedFluidPro:
O.engines +=
[PyRunner(command='activateFluidResol()',virtPeriod = 10,label = 'actFlResol')]
# else, it is a reloading of simulation, we load the necessary parameters and
the simulation
else :
execfile(scriptPath +'/paramSave.py') #Execute the file containing
the parameter registered previously
O.load(scriptPath + '/save.xml') #load the simulation where it
stopped
print '\n RELOAD THE PREVIOUS SIMULATION\n'
if imposedFluidPro==1 or initImposedFluidPro==1:
execfile('fluid.py') #Import the fluid profile (+turbStress if
necessary)
###################################################################################################################################
################################################### FUNCTION DEFINITION
#########################################################
###################################################################################################################################
######
######
### LET THE TIME FOR THE GRAVITY DEPOSITION AND ACTIVATE THE FLUID AT THE END
###
######
######
def gravityDeposition(lim):
#At the very start, unbalanced force can be low as there is only few
contacts, but it does not mean the packing is stable
#The rest will be run only if O.time>lim = estimation of the deposition
time
if O.time<lim : return
else :
newton.damping = 0.0
packing.dead = True # Remove the present engine for the
following
addFluidForce.dead = False # Activate the fluid force
fluidResol.dead = False # Activate the fluid profile resolution
# O.timingEnabled = True # Activate a function to check the time
taken by the different function
# If option activated, Activate the turbulent fluctuation
calculation
if turbModel==1:
turbFluct.dead = False
# If option activated, activate the function to perform
cumulated time averaging on DEM for fluidResol
if timeAveraging:
timeAve.dead = False
#If twoSize C++ option activated, express the size of the two
particles in the engine for later averaging
if twoSize and pythonAverage!=1:
addFluidForce.radiusPart1 = diameterPart1/2.
addFluidForce.radiusPart2 = diameterPart2/2.
addFluidForce.twoSize = True
if imposedFluidPro==1 or initImposedFluidPro==1:
fluidResol.dead = True #Force no fluid resolution
addFluidForce.vxFluid = vxFluid #pass the fluid
velocity vector to add drag force value
addFluidForce.simplifiedReynoldStresses = turbStress
#and reynolds stress tensor for fluctuations
else:
#Initialization of the fluid resolution
fluidResolution(fluidResolPeriod)
measurement.dead = False # Activate the Measurement
return
###############
#########################################
#############
###################
# Activate the fluid resolution. Used only when initImposedFluidPro option
activated ###
#############
###################
def activateFluidResol():
global initImposedFluidPro
if O.time>50:
fluidResol.dead = False #Activate the fluid resolution
initImposedFluidPro = 0 #Put the option to zero, to avoid
confusion in the script
actFlResol.dead = True #Kill the engine
#######################################################
###################################################################################################
############## Function to time-average the DEM data transmitted to fluid
resolution ###########
###################################################################################################
def timeAveragingDEMdata():
global phiPart_dt
global vxPart_dt
global averageDrag_dt
global phiPart1_dt
global averageDrag1_dt
global vxPart1_dt
global vyPart1_dt
global vzPart1_dt
global phiPart2_dt
global averageDrag2_dt
global vxPart2_dt
global vyPart2_dt
global vzPart2_dt
global n_dt
#Re-initialized at each new averaging step
if n_dt==0:
vxPart_dt = np.zeros(ndimz)
phiPart_dt = np.zeros(ndimz)
averageDrag_dt = np.zeros(ndimz)
if twoSize:
phiPart1_dt = np.zeros(ndimz)
averageDrag1_dt = np.zeros(ndimz)
vxPart1_dt = np.zeros(ndimz)
vyPart1_dt = np.zeros(ndimz)
vzPart1_dt = np.zeros(ndimz)
phiPart2_dt = np.zeros(ndimz)
averageDrag2_dt = np.zeros(ndimz)
vxPart2_dt = np.zeros(ndimz)
vyPart2_dt = np.zeros(ndimz)
vzPart2_dt = np.zeros(ndimz)
#Evaluate the average vx,vy,vz,phi,drag profiles and store it in
addFluidForce
addFluidForce.averageProfile()
#Store the values in
phiPart_dt += np.array(addFluidForce.phiPart)
vxPart_dt += np.array(addFluidForce.vxPart)
averageDrag_dt += np.array(addFluidForce.averageDrag)
if twoSize:
phiPart1_dt += np.array(addFluidForce.phiPart1)
averageDrag1_dt += np.array(addFluidForce.averageDrag1)
vxPart1_dt += np.array(addFluidForce.vxPart1)
vyPart1_dt += np.array(addFluidForce.vyPart1)
vzPart1_dt += np.array(addFluidForce.vzPart1)
phiPart2_dt += np.array(addFluidForce.phiPart2)
averageDrag2_dt += np.array(addFluidForce.averageDrag2)
vxPart2_dt += np.array(addFluidForce.vxPart2)
vyPart2_dt += np.array(addFluidForce.vyPart2)
vzPart2_dt += np.array(addFluidForce.vzPart2)
n_dt+=1.
###############
#########################################
################
#########################
########## EVALUATE THE TURBULENT FLUCTUATION ASSOCIATED TO EACH
PARTICLE##########
#################
#########################
# Argument: ndimz, turbulentViscosity, vxFluid, phiPart,
# Modification : turbFluct.virtPeriod, addFluidForce.bedElevation,
addFluidForce.simplifiedReynoldStresses, addFluidForce.velFluct
addFluidForce.simplifiedReynoldStresses = np.zeros(ndimz)
addFluidForce.bedElevation = 0.
addFluidForce.turbulentFluctuation()
addFluidForce.fluctTime = np.zeros(len(O.bodies))
def turbulentFluctuation():
global waterDepth
global turbStress
#For stability requirement at the initialization stage
if O.time<depoTime+0.5:
print 'No turbulent fluctuation in the initialization process
for stability reasons!'
turbFluct.virtPeriod = 0.5
else:
# Evaluate nBed, the position of the bed which is assumed to be
located around the first maximum of concentration when considering decreasing z.
nBed = 0.
bedElevation = 0.
for n in range(1,ndimz):
# if there is a peak and its value is superior to 0.5,
we consider it to be the position of the bed
if phiPart[ndimz - n] < phiPart[ndimz - n-1] and
phiPart[ndimz - n] > 0.5 :
nBed = ndimz - n
waterDepth = (ndimz-1 - nBed)*dz
bedElevation = fluidHeight - waterDepth
#Evaluate the bed elevation for the following
break
#(Re)Define the bed elevation over which fluid turbulent
fluctuations will be applied.
addFluidForce.bedElevation = bedElevation
#First turbulent fluctuation model: a unique constant lifetime
for the turbulent fluctuation, flucTimeScale
if turbModel1:
vMeanAboveBed = sum(vxFluid[nBed:])/(ndimz-nBed)
# fluid elocity scale in the water depth
flucTimeScale = waterDepth/vMeanAboveBed # time
scale of the fluctuation w_d/v, eddy turn over time
# New evaluation of the random fluid velocity
fluctuation for each particle.
addFluidForce.turbulentFluctuation()
turbFluct.virtPeriod = flucTimeScale #Actualize when
will be calculated the next fluctuations.
elif turbModel2:
# New evaluation of the random fluid velocity
fluctuation for each particle.
addFluidForce.turbulentFluctuationBIS()
turbFluct.virtPeriod =
min(addFluidForce.fluctTime)#Actualize when will be calculated the next fluct.
addFluidForce.dtFluct =
min(addFluidForce.fluctTime)#Actualize the time step of the function
if O.time<depoTime+1: #To avoid important virtPeriod at the
start, impose it for 1s. Fasten the transient.
turbFluct.virtPeriod = 0.1
if turbModel2: addFluidForce.dtFluct=0.1
###############
#########################################
###### ######
### CALL THE FLUID PROFILE RESOLUTION ROUTINE###
###### ######
#Arguments : ndimz, vxFluid, vxPart, dragTerm, densFluid, fluidHeight, sig,
dsig, diameterPart, phiPartFluid, kinematicViscoFluid, densPart, dpdx, slope,
tfin,d tFluid
#Modify : vxFluid, addFluidForce.vxFluid
testFluidResol = 0 # Test variable for the fluid resolution
def fluidResolution(tfin):
global phiPart
global vxPart
global vyPart
global vzPart
global vxFluid
global averageDrag
global turbulentViscosity
global testFluidResol
global dragTerm
global taufsi
global n_dt
global turbStress
if twoSize:
global averageDrag1
global phiPart1
global vxPart1
global vyPart1
global vzPart1
global averageDrag2
global phiPart2
global vxPart2
global vyPart2
global vzPart2
#Two different way of evaluating the average profiles: phiPart, vxPart,
averageDrag, (+ same for 1 and 2 with twoSize)
#From python
if pythonAverage==1 or debugStress==1:
averageProfile()
#Or from C++
elif timeAveraging:#with multiple time averaging made previously in
timeAveragingDEMdata
#Normalize the time averaged data
vxPart = vxPart_dt/n_dt
phiPart = phiPart_dt/n_dt
averageDrag = averageDrag_dt/n_dt
if twoSize:
phiPart1 = phiPart1_dt/n_dt
averageDrag1 = averageDrag1_dt/n_dt
vxPart1 = vxPart1_dt/n_dt
vyPart1 = vyPart1_dt/n_dt
vzPart1 = vzPart1_dt/n_dt
phiPart2 = phiPart2_dt/n_dt
averageDrag2 = averageDrag2_dt/n_dt
vxPart2 = vxPart2_dt/n_dt
vyPart2 = vyPart2_dt/n_dt
vzPart2 = vzPart2_dt/n_dt
#And report the other one (for usual measurement)
vyPart = np.array(addFluidForce.vyPart)
vzPart = np.array(addFluidForce.vzPart)
n_dt = 0
else:#with simple time averaging
if raphaelComputer!=1:
#Evaluate the average vx,vy,vz,phi,drag profiles and
store it in addFluidForce
addFluidForce.averageProfile()
else:#On old ubuntu/Yade version, needs to do it differently:
#At next time step, Evaluate the average
vx,vy,vz,phi,drag profiles and store it in addFluidForce
if testFluidResol == 0:
addFluidForce.activateAverage = True
#Activate the average for next time step
fluidResol.iterPeriod = 1 #For it to be
re-executed at next time step
testFluidResol = 1
return
## Import the necessary average variables for the fluid
resolution from addFluidForce : phiPart,vxPart, averageDrag
vxPart = np.array(addFluidForce.vxPart)
vyPart = np.array(addFluidForce.vyPart)
vzPart = np.array(addFluidForce.vzPart)
phiPart = np.array(addFluidForce.phiPart)
averageDrag = np.array(addFluidForce.averageDrag)
#Back to the usual resolution period for next step
fluidResol.iterPeriod = -1 #For it not to be re-executed
next time
fluidResol.virtPeriod = fluidResolPeriod #Back to the
normal period
testFluidResol = 0 #Re-initialize the test
#If two particle size in the simulation, store complementary
necessary average
if twoSize:
vCell = dz*width*length
#C++ function evaluating the depth profiles of solid
volume fraction and solid velocities
phiPart1=np.array(addFluidForce.phiPart1)
phiPart2=np.array(addFluidForce.phiPart2)
averageDrag1=np.array(addFluidForce.averageDrag1)
averageDrag2=np.array(addFluidForce.averageDrag2)
vxPart1=np.array(addFluidForce.vxPart1)
vxPart2=np.array(addFluidForce.vxPart2)
vyPart1=np.array(addFluidForce.vyPart1)
vyPart2=np.array(addFluidForce.vyPart2)
vzPart1=np.array(addFluidForce.vzPart1)
vzPart2=np.array(addFluidForce.vzPart2)
#
## Prepare the variables for the fluid resolution
#
# Threshold the solid vol. fraction in the bed to remove the
oscillations due to layering and avoid turbulence in the bed
phiPartFluid = np.zeros(ndimz)
for n in range(1,ndimz+1):
# If the concentration inside the bed reach the maximum
concentration (imposed in the fluid solver). Then all the values of
concentration under this one are going to be thresholded at phiPartMax, in
order to damp the turbulence inside the bed.
if phiPart[ndimz - n] > phiPartMax:
for i in range(n,ndimz+1):
phiPartFluid[ndimz - i] = phiPartMax
break
#Otherwise we report the normal value.
phiPartFluid[ndimz - n] = phiPart[ndimz - n]
#Density of particle per unit volume N/Vcell:
nPart = phiPart[:ndimz]/partVolume
#Calculate the momentum transfer from the drag force interaction :
N<fd>/Vcell
dragTerm = nPart*averageDrag[:ndimz]
if twoSize: #if two size, sum the two contributions n1<f1> + n2 <f2>,
formulation necessary to keep energy conservation.
dragTerm = phiPart1[:ndimz]/partVolume1*averageDrag1[:ndimz] +
phiPart2[:ndimz]/partVolume2*averageDrag2[:ndimz]
#Create Taufsi/rhof = dragTerm/(rhof(vf-vxp)) to transmit to the fluid
code
taufsi = np.zeros(ndimz)
for i in range(1,ndimz):
urel = max(abs(vxFluid[i] - vxPart[i]),1e-5) #Limit the
value to avoid division by 0
#Limit the max value of taufsi to the fluid resolution limit,
i.e. 1/(10dt) and required positive (charact. time)
taufsi[i] = max(0,min(dragTerm[i]/urel/densFluid,
pow(10*dtFluid,-1)))
#
## Call the routine for fluid resolution
#
[vxFluid_np,turbulentViscosity] =
nsmp1d_yade(fluidHeight,sig,dsig,diameterPart,vxFluid,1e0-phiPartFluid[0:ndimz],densFluid,kinematicViscoFluid,vxPart[0:ndimz],phiPartFluid[0:ndimz],densPart,dpdx,slope,tfin,dtFluid,taufsi)
#Debugstress option
if debugStress==1:
global taufsiBIS
taufsiBIS = np.zeros(ndimz)
for i in range(1,ndimz):
urel = max(abs(vxFluid[i] - vxPart[i]),1e-5)
taufsiBIS[i] = min(dragTermBIS[i]/urel/densFluid,
pow(10*dtFluid,-1))
global vxFluid_before
vxFluid_before = vxFluid
vxFluid = vxFluid_np
path = scriptPath + '/'+ str(fileNumber)+'_debug.py'
globalParam =
['vxPart','vxPartBIS','phiPart','phiPartBIS','vxFluid','vxFluid_before',
'averageDrag','averageDragBIS','dragTerm','dragTermBIS','taufsi','taufsiBIS']
Save(path, globalParam)
#update
vxFluid = vxFluid_np
#Modify the velocity in the HydroForceEngine which apply the force
addFluidForce.vxFluid = vxFluid
#Evaluate the turbulent shear stress for the turbulent fluctuations
model and Shields number evaluation
turbStress = np.zeros(ndimz)
n = 0
for nu_t in turbulentViscosity[:-2]: # Turbulent viscosity is
global, calculated in fluidResolution
turbStress[n] = nu_t*(vxFluid[n+1]-vxFluid[n])/dz #tau/rho
n+=1
addFluidForce.simplifiedReynoldStresses = turbStress #Actualize the
Reynolds stress tensor (/rho)value
###############
#########################################
####### ########
### OUTPUT ###
####### ########
#Function which register the different quantities asked for post processing
#Arguments : len(O.bodies), O.bodies, scriptPath, fileNumber, measurePeriod,
savePeriod, qsMean, phiPart, vxPart, vxFluid, turbulentViscosity, dragTerm,
waterDepth
#Modify : fileNumber, 'paramSave.py', 'save.xml'
def measure():
global fileNumber
global qsMean
global stressProfile
global vxPart
global vyPart
global vzPart
global phiPart
global averageDrag
global phiPart1
global phiPart2
global averageDrag1
global averageDrag2
global vxPart1
global vxPart2
global vyPart1
global vyPart2
global vzPart1
global vzPart2
if imposedFluidPro==1 or initImposedFluidPro==1:#Then you should
evaluate the average (as not done in fluidResolution)
if pythonAverage==1 or debugStress==1:
averageProfile()
#Or from C++
else:
#Evaluate the average vx,vy,vz,phi,drag profiles and
store it in addFluidForce
addFluidForce.averageProfile()
#Replace it in the global python variables
vxPart = np.array(addFluidForce.vxPart)
vyPart = np.array(addFluidForce.vyPart)
vzPart = np.array(addFluidForce.vzPart)
phiPart = np.array(addFluidForce.phiPart)
averageDrag = np.array(addFluidForce.averageDrag)
#If two particle size in the simulation, store
complementary necessary average
if twoSize:
vCell = dz*width*length
#C++ function evaluating the depth profiles of
solid volume fraction and solid velocities
phiPart1=np.array(addFluidForce.phiPart1)
phiPart2=np.array(addFluidForce.phiPart2)
averageDrag1=np.array(addFluidForce.averageDrag1)
averageDrag2=np.array(addFluidForce.averageDrag2)
vxPart1=np.array(addFluidForce.vxPart1)
vxPart2=np.array(addFluidForce.vxPart2)
vyPart1=np.array(addFluidForce.vyPart1)
vyPart2=np.array(addFluidForce.vyPart2)
vzPart1=np.array(addFluidForce.vzPart1)
vzPart2=np.array(addFluidForce.vzPart2)
# Save all the following variables/split in two file for faster post
processing
if stressSave == 1:
if pythonStressCalculation:
stressProfile =
stressFieldCalculation(ndimz,dz*width*length,vxPart,dz)
else:
stressProfile =
getStressProfile(dz*width*length,ndimz,dz,groundPosition,vxPart,vyPart,vzPart)
for n in range(0,ndimz):
stressProfile[0][n] =
np.array(stressProfile[0][n]).reshape(1,9)[0]
stressProfile[1][n] =
np.array(stressProfile[1][n]).reshape(1,9)[0]
globalParam = ['stressProfile']
path = scriptPath + '/'+ str(fileNumber)+'_stress.py'
Save(path, globalParam)
#if necessary output file with the positions and velocities of the
particles
if 0:
global positionPartX
global positionPartZ
global velocityPartX
global velocityPartZ
positionPartX = np.zeros(len(O.bodies)) #Defined global in
order to use
positionPartZ = np.zeros(len(O.bodies)) #Defined global in
order to use
velocityPartX = np.zeros(len(O.bodies)) # the saving function
velocityPartZ = np.zeros(len(O.bodies)) # the saving function
for b in O.bodies:
if isinstance(b.shape,Sphere): # If it is a sphere
positionPartX[b.id] = b.state.pos[0]
positionPartZ[b.id] = b.state.pos[2]
velocityPartX[b.id] = b.state.vel[0]
velocityPartZ[b.id] = b.state.vel[2]
path = scriptPath + '/'+ str(fileNumber)+'_2.py'
globalParam = ['positionPartX',
'positionPartZ','velocityPartX', 'velocityPartZ']
Save(path, globalParam)
#Evaluate the dimensionless sediment rate to save
qsMean = sum(phiPart*vxPart)*dz/sqrt((densPart/densFluid -
1)*abs(gravityVector[2])*pow(diameterPart,3))
path = scriptPath + '/'+ str(fileNumber)+'.py'
globalParam =
['qsMean','phiPart','vxPart','vyPart','vzPart','vxFluid','turbulentViscosity','dragTerm','turbStress']
if twoSize:
globalParam
=['qsMean','phiPart','vxPart','vyPart','vzPart','vxFluid','turbulentViscosity','dragTerm','turbStress',
'phiPart1','vxPart1','vyPart1','vzPart1','phiPart2','vxPart2','vyPart2','vzPart2']
Save(path, globalParam)
fileNumber +=1
#PHF if twoSize!=1 and O.time>=150:
# if O.time>=150:
# print '\n End of the simulation, simulated 150s as required !\n
'
# O.pause()
# O.exitNoBacktrace()
#For reload, save sometimes the variables necessary to reload. (does
not work anymore...)
if (fileNumber*measurePeriod)%savePeriod == 0:
#Save the simulation of yade the global python parameter in
order to be able to recover the simulation
O.save(scriptPath + '/'+'save.xml')
# Save the global parameters which shouldn't be initialized to
zero at the reloading.
globalParam =
['phiPart','fileNumber','vxFluid','turbulentViscosity','waterDepth','turbStress']
Save(scriptPath + '/'+'paramSave.py', globalParam)
##
# Evaluation and printing of the Shields number or qsMean.
if imposedFluidPro==1 or initImposedFluidPro==1:
print 'qsMean', qsMean
else:
if twoSize:
#Consider the average diameter in the layer in motion
(which is the one to consider for Shields number)
dummy = np.where(phiPart[:ndimz]<0.8*phiPartMax)
phi1Sum = sum(phiPart1[dummy])
phi2Sum = sum(phiPart2[dummy])
#Volume weighted average diameter: N1Vp1 d1 + N2Vp2
d2/(N1Vp1 + N2Vp2) (*VCell/VCell)
diameterMix = (diameterPart1*phi1Sum +
diameterPart2*phi2Sum)/(phi1Sum + phi2Sum)
else:
#Only one size in the simulation...
diameterMix = diameterPart
shieldsNumber =
densFluid*max(turbStress)/((densPart-densFluid)*diameterMix*abs(gravityVector[2]))
print 'Shields number', shieldsNumber
###############
#########################################
##########
##########
### EVALUATE THE MEAN VOLUME FRACTION, DRAG AND STREAMWISE VELOCITY OF THE
PARTICLES IN FUNCTION OF THE DEPTH###
##########
##########
#Arguments : ndimz, dz, length, Vlayer, groundPosition, radius, phiPart
#Modify : phiPart, vxPart, averageDrag
def averageProfile():
global phiPart
global averageDrag
global vxPart
global vyPart
global vzPart
global qsMean
global phiPart1
global averageDrag1
global phiPart2
global averageDrag2
#Initialization of the concentration and the streamwise velocity. (two
times the water height in order to take into account the particles which are
also above the fluid in the concentration and velocity profile
phiPart_loc = np.zeros(ndimz) #Mean particle concentration in
function of depth
vxPart_loc = np.zeros(ndimz) #Mean particle velocity in
function of depth
vyPart_loc = np.zeros(ndimz) #Mean particle velocity in
function of depth
vzPart_loc = np.zeros(ndimz) #Mean particle velocity in
function of depth
averageDrag_loc = np.zeros(ndimz) #Mean drag
phiPart1_loc = np.zeros(ndimz)
averageDrag1_loc = np.zeros(ndimz)
phiPart2_loc = np.zeros(ndimz)
averageDrag2_loc = np.zeros(ndimz)
if debugStress==1:
global vxPartBIS
global phiPartBIS
global averageDragBIS
global dragTermBIS
averageDragBIS_loc = np.zeros(ndimz)
phiPartBIS_loc = np.zeros(ndimz)
vxPartBIS_loc = np.zeros(ndimz)
dragTermBIS_loc = np.zeros(ndimz)
#Import the vector of velocity fluctuation associated to each particle
in C++ code
VxFluct = addFluidForce.vFluctX
VzFluct = addFluidForce.vFluctZ
shiftPosition = 1.2*diameterPart #To shift the position
in order to avoid to have negative position (not compatible with the
volume calculation used later)
## Loop over the file to determine the concentration and the mean
velocity
for b in O.bodies:
if isinstance(b.shape,Sphere): # If it is a sphere
z = b.state.pos[2]-groundPosition # Altitude of
the particle redefined with the position of the bottom wall as the
radius = b.shape.radius
Np = int(z/dz) #Layer corresponding to
the position of the center of the particle
#Evaluate the drag force undergone by the particle
(Dallavalle formulation + Richardson Zaki correction)
uRel = Vector3(0.0,.0,0.0)
if Np<ndimz and Np>=0 : uRel =
Vector3(vxFluid[Np]+VxFluct[b.id],0.0,VzFluct[b.id]) - b.state.vel
fDrag =
0.5*pi*pow(radius,2)*densFluid*(0.44*uRel.norm()+24.4*kinematicViscoFluid/(2*radius))*pow(1-phiPart[Np],-expoRZ)*uRel
if debugStress==1:
dragTermBIS_loc[Np]+=fDrag[0]
phiPartBIS_loc[Np]+=1.0
vxPartBIS_loc[Np]+=b.state.vel[0]
#Shift the position to avoid z<0 which induces some
problems in the volume calculation. This is done here in order not to change
the value of the drag. Evaluate all the parameters after this shift. (Np in
particular is reevaluated
z+=shiftPosition
Np = int(z/dz) #new evaluation of Np
after the shift. DO NOT REMOVE, IMPORTANT
minZ = int((z-radius)/dz) #Lowest layer occupied
by the particle
maxZ = int((z+radius)/dz) #Highest layer occupied
by the particle
deltaCenter = z - Np*dz #length between the
real center of the particle and the position of the layer in which it is
considered to be (Np).
#loop from the bottom to the top layer inside which the
particle is contained
numLayer = minZ
while numLayer <= maxZ:
if numLayer>=0 and numLayer<ndimz:
zInf = (numLayer-Np-1)*dz + deltaCenter
zSup = (numLayer-Np)*dz + deltaCenter
if numLayer==minZ: zInf = -radius
if numLayer==maxZ: zSup = radius
#evaluate the volume of the part of the
particle contained in the layer considered (the formula results from the
analytical resolution of the integral of a layer of sphere in cylindrical
coordinate)
V = math.pi*radius*radius*(zSup-zInf
+(pow(zInf,3)-pow(zSup,3))/(3*radius*radius))
#Add the quantities to the (not
normalized for now) mean particle velocity and concentration
phiPart_loc[numLayer]+=V
vxPart_loc[numLayer]+=V*b.state.vel[0]
vyPart_loc[numLayer]+=V*b.state.vel[1]
vzPart_loc[numLayer]+=V*b.state.vel[2]
averageDrag_loc[numLayer]+=V*fDrag[0]
if twoSize:
if radius==diameterPart1/2.0:
phiPart1_loc[numLayer]+=V
averageDrag1_loc[numLayer]+=V*fDrag[0]
elif radius==diameterPart2/2.0:
phiPart2_loc[numLayer]+=V
averageDrag2_loc[numLayer]+=V*fDrag[0]
#Increment for the next loop
numLayer+=1
#Normalize the weighted velocity/different mean quantities of each
layer by the total particle volume in this layer
n = 0
for i in phiPart_loc:
if i: #if non zero, to avoid division by zero
vxPart_loc[n]/=i
vyPart_loc[n]/=i
vzPart_loc[n]/=i
averageDrag_loc[n]/=i
if twoSize:
if phiPart1_loc[n]:
averageDrag1_loc[n]/=phiPart1_loc[n]
if phiPart2_loc[n]:
averageDrag2_loc[n]/=phiPart2_loc[n]
n += 1
#Normalization of the concentration
phiPart_loc/=VLayer
phiPart1_loc/=VLayer
phiPart2_loc/=VLayer
if debugStress==1:
#Normalization
n = 0
for i in phiPartBIS_loc:#Loop over the number of part contained
in each layer
#Normalize the average drag and particle velocity
if i:
averageDragBIS_loc[n] = dragTermBIS_loc[n]/i
vxPartBIS_loc[n]/=i
n+=1
phiPartBIS_loc*=(partVolume/VLayer)
dragTermBIS_loc/=VLayer
#Report the value in global variables
vxPartBIS = vxPartBIS_loc
phiPartBIS = phiPartBIS_loc
averageDragBIS = averageDragBIS_loc
dragTermBIS = dragTermBIS_loc
#Resize the vector to remove the effect of the shift previously
introduced
phiPart_loc = Resize(phiPart_loc,int(shiftPosition/dz))
vxPart_loc = Resize(vxPart_loc,int(shiftPosition/dz))
vyPart_loc = Resize(vyPart_loc,int(shiftPosition/dz))
vzPart_loc = Resize(vzPart_loc,int(shiftPosition/dz))
averageDrag_loc = Resize(averageDrag_loc,int(shiftPosition/dz))
phiPart1_loc = Resize(phiPart1_loc,int(shiftPosition/dz))
averageDrag1_loc = Resize(averageDrag1_loc,int(shiftPosition/dz))
phiPart2_loc = Resize(phiPart2_loc,int(shiftPosition/dz))
averageDrag2_loc = Resize(averageDrag2_loc,int(shiftPosition/dz))
#Impose variables at the bottom to be consistent with the fluid
formulation
averageDrag_loc[0] = 0.0
vxPart_loc[0] = 0.0
###Actualize the global variables
phiPart = phiPart_loc
vxPart = vxPart_loc
vyPart = vyPart_loc
vzPart = vzPart_loc
averageDrag = averageDrag_loc
phiPart1 = phiPart1_loc
averageDrag1 = averageDrag1_loc
phiPart2 = phiPart2_loc
averageDrag2 = averageDrag2_loc
#And the one associated to the fluid force application
addFluidForce.phiPart = phiPart_loc
return
###############
#########################################
###### ######
## Function to resize the vector removing the n first case##
###### ######
def Resize(a,n):
import numpy as np
l = len(a)
b = np.zeros(l)
k = 0
for i in range(n,l):
b[k] = a[i]
k+=1
return b
###############
#########################################
###### ######
### Save global variables ###
###### ######
#Function to save global variables in a python file which can be re-executed
for later
def Save(filePathName, globalVariables):
f = open(filePathName,'w')
f.write('from numpy import *\n')
for i in globalVariables:
f.write(i + ' = '+repr(globals()[i]) + '\n')
f.close()
###############
#########################################
###### ######
### python Stress field calculation ###
###### ######
def stressFieldCalculation(ndimz,vCell,vxPart,dz):
#Initialization
granularTemp_loc = np.zeros(ndimz)
N_loc = np.zeros(ndimz)
stressTensorField_loc = []
kineticStressTensorField_loc = []
for i in range(0,ndimz):
stressTensorField_loc += [Matrix3(0,0,0,0,0,0,0,0,0)]
kineticStressTensorField_loc += [Matrix3(0,0,0,0,0,0,0,0,0)]
###
## Dynamical stress tensor part
###
for b in O.bodies:
Np = int((b.state.pos[2] - groundPosition)/dz)
if Np>=0 and Np<ndimz:
vFluct = b.state.vel -
Vector3(vxPart[Np],vyPart[Np],vzPart[Np])
stressTensorField_loc[Np]+=
-1/vCell*b.state.mass*vFluct.outer(vFluct)
kineticStressTensorField_loc[Np]+=
-1/vCell*b.state.mass*vFluct.outer(vFluct)
granularTemp_loc[Np] += 1/3.0*(pow(vFluct[0],2) +
pow(vFluct[1],2) + pow(vFluct[2],2))
N_loc[Np]+=1.0
###
## Love-Weber stress tensor part (contact
###
for Int in O.interactions:
if O.bodies[Int.id1].dynamic or O.bodies[Int.id2].dynamic:
Np1 = int((O.bodies[Int.id1].state.pos[2] -
groundPosition)/dz)
Np2 = int((O.bodies[Int.id2].state.pos[2] -
groundPosition)/dz)
#Vector joining the two particles (from id2 to id1)
x12 = O.bodies[Int.id1].state.pos -
O.bodies[Int.id2].state.pos
#Contact vector of the two particles (defined as
applied by id1 on id2)
fContact = Int.phys.normalForce + Int.phys.shearForce
#Remove the artificial effect of the periodic BC if one
particle goes to the other side
x12 -= O.cell.hSize*Int.cellDist
#No reason to have one upon another, so check and
define the min and max correpondingly
if Np2==Np1: #If the contact vector is entirely
inside the layer, add it directly.
if Np1>=0 and Np1<ndimz:
stressTensorField_loc[Np1]+=
1/vCell*fContact.outer(x12) #outer product of contact force and branch
vector
else:
if Np1>Np2: #Particle 1 above 2 (z
component)
minZ = Np2
minPosZ =
O.bodies[Int.id2].state.pos[2] - groundPosition
maxZ = Np1
maxPosZ =
O.bodies[Int.id1].state.pos[2] - groundPosition
elif Np2>Np1: #Particle 2 above 1 (z
component)
minZ = Np1
minPosZ =
O.bodies[Int.id1].state.pos[2] - groundPosition
maxZ = Np2
maxPosZ =
O.bodies[Int.id2].state.pos[2] - groundPosition
#Loop over the layers containing the vector
connecting the two particles
numLayer = minZ
while numLayer<=maxZ:
if numLayer>=0 and numLayer<ndimz:
deltaZ = dz
if numLayer==minZ:
deltaZ = dz - (minPosZ -
minZ*dz)
elif numLayer==maxZ:
deltaZ = maxPosZ -
maxZ*dz
branchVectCell=deltaZ/abs(x12[2])*Vector3(x12[0],x12[1],x12[2])
#Add to stress tensor of the
considered layer the part considered
stressTensorField_loc[numLayer]+= 1/vCell*fContact.outer(branchVectCell)
#outer product
#Increment the layer
numLayer+=1
#Normalization
for n in range(0,len(N_loc)):
if N_loc[n]>0:
granularTemp_loc[n]/=N_loc[n]
return
[stressTensorField_loc,kineticStressTensorField_loc,granularTemp_loc]
###############
#########################################
###DEBUG/COMPARISON WITH getStress
#To obtain exactly the same as getStress, you need to shift groundPosition ->
groundPosition - diameterPart + to remove the condition on the dynamics of the
particles, to avoid the non contribution of the bottom wall and particles.
#This function might be also useful, it is reproducing exactly getStress in
python
def getStressBIS(vCell):
stressTensorField = Matrix3(0,0,0,0,0,0,0,0,0)
for Int in O.interactions:
x12 = O.bodies[Int.id1].state.pos - O.bodies[Int.id2].state.pos
x12 -= O.cell.hSize*Int.cellDist
fContact = Int.phys.normalForce + Int.phys.shearForce
stressTensorField+= fContact.outer(x12) #outer product of
contact force a
return stressTensorField/vCell
###############
#########################################
--
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
