yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3919
committer: Raphael Maurin <raph_maurin@xxxxxxxxxxx>
timestamp: Sat 2016-07-09 00:21:49 +0200
message:
  Add two example scripts for the use of HydroForceEngine
  One simple example to get familiar with the application of buoyancy (buoyantParticles.py).
  One example simulating a fluidized bed configuration to get familiar with both the application of a fluid velocity profile and of turbulent fluctuation with the DRW models (fluidizedBed.py)
added:
  examples/HydroForceEngine/buoyantParticles.py
  examples/HydroForceEngine/fluidizedBed.py


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== added file 'examples/HydroForceEngine/buoyantParticles.py'
--- examples/HydroForceEngine/buoyantParticles.py	1970-01-01 00:00:00 +0000
+++ examples/HydroForceEngine/buoyantParticles.py	2016-07-08 22:21:49 +0000
@@ -0,0 +1,104 @@
+#########################################################################################################################################################################
+# Author: Raphael Maurin, raphael.maurin@xxxxxxx
+# 07/07/2016
+#
+# Example script to use HydroForceEngine in order to apply buoyancy force to particles. 
+# The fluid is supposed at rest (vxFluidPY is a zero vector) so that the particles are only submitted to a drag force opposing the motion of the particle. 
+# Three spheres with density 1500, 1000 and 500kg/m3 are positionned at the middle of a fluid sample of density 1000 kg/m3, and let evolved with time
+# We observe clearly that one sphere sediment down to the bottom, another one is rising to the top of the water free-surface and the third one does not move
+#
+############################################################################################################################################################################
+
+#Import libraries
+from yade import pack, plot
+import math
+import random as rand
+import numpy as np
+
+
+##
+## Main parameters of the simulation
+##
+
+#Particles
+diameterPart = 6e-3	#Diameter of the particles, in meter
+restitCoef = 0.5	#Restitution coefficient of the particles
+partFrictAngle = atan(0.4)	#friction angle of the particles, in radian
+densPart1 = 1050		#density of the particles, in kg/m3
+densPart2 = 1000		#density of the particles, in kg/m3
+densPart3 = 950		#density of the particles, in kg/m3
+
+#fluid
+densFluidPY = 1000.	#Density of the fluid
+kinematicViscoFluid = 1e-6	#kinematic viscosity of the fluid
+fluidHeight = 20.*diameterPart	#Water depth in m.
+ndimz = 1
+groundPosition = 0.#Definition of the position of the ground, in m
+gravityVector = Vector3(0,0.0,-9.81) #Gravity vector
+
+#Particles contact law/material parameters
+normalStiffness = 1e4
+youngMod = normalStiffness/diameterPart	#Young modulus of the particles from the stiffness wanted.
+poissonRatio = 0.5	#poisson's ratio of the particles. Classical values, does not have much influence
+
+#Material of particle 1, 2, and 3, with different density densPart1, densPart2 and densPart3 defined above (1500, 1000 and 500kg/m3 by default) 
+O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod, poisson=poissonRatio, density=densPart1, frictionAngle=partFrictAngle, label='Mat1'))  
+O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod, poisson=poissonRatio, density=densPart2, frictionAngle=partFrictAngle, label='Mat2'))  
+O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod, poisson=poissonRatio, density=densPart3, frictionAngle=partFrictAngle, label='Mat3'))  
+
+#Time of simulation
+endTime = 2.
+
+########################
+## FRAMEWORK CREATION ##
+########################
+
+# Reference walls: build a wall at the ground and draw the position of the free-surface to have a reference for the eyes in the 3D view
+lowPlane = box(center= (0, 0,groundPosition),extents=(200,200,0),fixed=True,wire=False,color = (0.,1.,0.))
+WaterSurface = box(center= (0,0,groundPosition+fluidHeight),extents=(200,200,0),fixed=True,wire=False,color = (0,0,1),mask = 0)
+O.bodies.append([lowPlane,WaterSurface]) #add to simulation
+
+id1 = O.bodies.append(sphere(center = (0,2*diameterPart,groundPosition + fluidHeight/2.), radius = diameterPart/2.,material = 'Mat1'))
+id2 = O.bodies.append(sphere(center = (0,0,groundPosition + fluidHeight/2.), radius = diameterPart/2.,material = 'Mat2'))
+id3 = O.bodies.append(sphere(center = (0,-2*diameterPart,groundPosition + fluidHeight/2.), radius = diameterPart/2.,material = 'Mat3'))
+
+# Collect the ids of the spheres which are dynamic to add a fluid force through HydroForceEngines
+idApplyForce = [id1,id2,id3]
+
+
+#########################
+#### 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 = densFluidPY,viscoDyn = kinematicViscoFluid*densFluidPY,zRef = groundPosition,gravity = gravityVector,deltaZ = fluidHeight/ndimz,expoRZ = 0.,lift = False,nCell = ndimz,vCell = 1.,vxFluid = np.zeros(ndimz),phiPart = np.zeros(ndimz),vxPart = np.zeros(ndimz),vFluctX = np.zeros(len(O.bodies)),vFluctY = np.zeros(len(O.bodies)),vFluctZ = np.zeros(len(O.bodies)),ids = idApplyForce, label = 'hydroEngine'),
+	#To plot the wall normal position of the spheres with time
+	PyRunner(command = 'plotPos()', virtPeriod = 0.01, label = 'plot'),
+	# Integrate the equation and calculate the new position/velocities...
+	NewtonIntegrator(gravity=gravityVector, label='newtonIntegr')
+	]
+#save the initial configuration to be able to recharge the simulation starting configuration easily
+O.saveTmp()
+#Time step
+O.dt = 5e-7 #Low no purpose, in order to observe the sedimentation
+
+#Plot the wall normal position of the spheres with time
+def plotPos():
+	plot.addData(z1 = O.bodies[2].state.pos[2]/fluidHeight,z2 = O.bodies[3].state.pos[2]/fluidHeight,z3 = O.bodies[4].state.pos[2]/fluidHeight, time = O.time)
+	if O.time>endTime:
+		print('\nEnd of the simulation, {0}s simulated as asked!\n'.format(endTime))
+		O.pause()
+
+plot.plots={'time':('z1','z2','z3')}
+plot.plot()

=== added file 'examples/HydroForceEngine/fluidizedBed.py'
--- examples/HydroForceEngine/fluidizedBed.py	1970-01-01 00:00:00 +0000
+++ examples/HydroForceEngine/fluidizedBed.py	2016-07-08 22:21:49 +0000
@@ -0,0 +1,172 @@
+#########################################################################################################################################################################
+# Author: Raphael Maurin, raphael.maurin@xxxxxxx
+# 08/07/2016
+#
+#  Very simplified fluidized bed simulations in order to underline the possibility of HydroForceEngine. 
+#  Particles are deposited under gravity inside a box. Once the particles at rest, a constant fluid velocity is applied against gravity, submitting the
+#  particles to a constant drag force. Then, a discrete random walk model is applied to mimick the effect of the turbulent fluid velocity fluctuations. 
+#  It associates to each particle a random fluid velocity fluctuations in the x, y and z directions, which are taken into account in the evaluation
+#  of the drag applied by the fluid on the particle. The intensity of the fluctuation is based on the value of u*, which is imposed through simplifiedReynoldsStress = u*^2
+#  The values taken for the fluid velocity and simplifiedReynolds stress have been arbitrarily chosen to have a nice rendering
+#
+#  The example allows to get familiar with HydroForceEngine and in particular with the included DRW model functions (turbulentFluctuationFluidizedBed() here)
+#  For details on the process, it is necessary to have a look to the documentation and the C++ code in pkg/common/ForceEngine.cpp and hpp as HydroForceEngine contains
+#  multiple parameters
+############################################################################################################################################################################
+
+# ATTENTION: to fit the formulation of HydroForceEngine, the fluid velocity can only be applied along the x axis. Therefore, the gravity is here aligned with x, and all the 
+# configuration is defined accordingly. For the 3D visualization, two walls of the cell have been made transparent and clicking on the right xyz button after openning the 
+# 3D view allows one to see the sample in the usual way with gravity going from top to bottom. 
+
+
+#Import libraries
+from yade import pack, plot
+import math
+import random as rand
+import numpy as np
+
+
+##
+## Main parameters of the simulation
+##
+
+#Particles
+diameterPart = 6e-3	#Diameter of the particles, in meter
+densPart = 2500		#density of the particles, in kg/m3
+restitCoef = 0.8	#Restitution coefficient of the particles
+partFrictAngle = atan(0.4)	#friction angle of the particles, in radian
+
+#fluid
+densFluidPY = 1.225	#Density of the fluid, (air) in kg/m3
+kinematicViscoFluid = 1.48e-5	#kinematic viscosity of the fluid, (air) in m2/s
+
+#Configuration: inclined channel
+lengthCell = 100	#length cell along the gravity axis (x), in diameter
+widthCell = 10	#length cell along the two other axis, in diameter
+Nlayer = 2.5	#nb of layer of particle deposited, in diameter
+phiPartMax = 0.61	#Value of the dense packing solid volume fraction
+
+endTime = 10	#Time simulated (in seconds)
+
+
+##
+## Secondary parameters of the simulation
+##
+
+expoDrag_PY = 3.1	# Richardson Zaki exponent for the hindrance function of the drag force applied to the particles
+ndimz = 20	#Number of cells in the height
+dz =  widthCell*diameterPart/ndimz	# Fluid discretization step in the wall-normal direction	
+
+# Initialization of the main vectors
+vxFluidPY = np.ones(ndimz)*18.5	# Vertical fluid velocity profile: u^f = u_x^f(z) e_x, with x the streamwise direction and z the wall-normal
+
+phiPartPY = np.zeros(ndimz)	# Vertical particle volume fraction profile
+vxPartPY = np.zeros(ndimz)	# Vertical average particle velocity profile
+
+#Geometrical configuration, define useful quantities
+length = lengthCell*diameterPart #length of the stream, in m
+width = widthCell*diameterPart  #width of the stream, in m
+groundPosition = 0. #Definition of the position of the ground, in m
+gravityVector = Vector3(-9.81,0.0,0.) #Gravity vector. Inclined along x to have an effect opposed to the fluid flow (which can only be along x)
+
+#Particles contact law/material parameters
+maxPressure = (densPart-densFluidPY)*phiPartMax*Nlayer*diameterPart*9.81#Estimated max particle pressure from the static load
+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.
+poissonRatio = 0.5	#poisson's ratio of the particles. Classical values, does not have much influence
+O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod, poisson=poissonRatio, density=densPart, frictionAngle=partFrictAngle, label='Mat'))  
+
+
+########################
+## FRAMEWORK CREATION ##
+########################
+
+# Walls to create a box
+lowPlane = box(center= (groundPosition, width/2.0,width/2.),extents=(0,width/2.,width/2.),fixed=True,wire=False,color = (0.,1.,0.),material = 'Mat')
+sidePlane1 = box(center= (length/2., width,width/2.),extents=(length/2.,0.,width/2.),fixed=True,wire=True,color = (0.,1.,0.),material = 'Mat')
+sidePlane2 = box(center= (length/2., 0.,width/2.),extents=(length/2.,0.,width/2.),fixed=True,wire=True,color = (0.,1.,0.),material = 'Mat')
+sidePlane3 = box(center= (length/2.,width/2., 0.),extents=(length/2.,width/2.,0.),fixed=True,wire=False,color = (0.,1.,0.),material = 'Mat')	#Made invisible (wire = True) in order to see inside the cell
+sidePlane4 = box(center= (length/2.,width/2., width),extents=(length/2.,width/2.,0.),fixed=True,wire=False,color = (0.,1.,0.),material = 'Mat')	#Made invisible (wire = True) in order to see inside the cell
+O.bodies.append([lowPlane,sidePlane1,sidePlane2,sidePlane3,sidePlane4])
+
+#Create a loose cloud of particle inside the cell
+partCloud = pack.SpherePack()
+partVolume = pi/6.*pow(diameterPart,3) #Volume of a particle
+partNumber = int(Nlayer*phiPartMax*diameterPart*width*width/partVolume) #Volume of beads to obtain Nlayer layers of particles
+partCloud.makeCloud(minCorner=(0,0.,0),maxCorner=(length,width,width),rRelFuzz=0., rMean=diameterPart/2.0, num = partNumber)
+partCloud.toSimulation(material='Mat') #Send this packing to simulation with material Mat
+#Evaluate the deposition time considering the free-fall time of the highest particle to the ground
+depoTime = sqrt(length*2/9.31)
+
+# Collect the ids of the spheres which are dynamic to add a fluid force through HydroForceEngines
+idApplyForce = []
+for b in O.bodies: 
+	if isinstance(b.shape,Sphere) and b.dynamic:
+		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 = densFluidPY,viscoDyn = kinematicViscoFluid*densFluidPY,zRef = groundPosition,gravity = gravityVector,deltaZ = dz,expoRZ = expoDrag_PY,lift = False,nCell = ndimz,vCell = length*width*dz ,vxFluid = vxFluidPY,phiPart = phiPartPY,vxPart = vxPartPY,ids = idApplyForce,vFluctX = np.zeros(len(O.bodies)),vFluctY = np.zeros(len(O.bodies)),vFluctZ = np.zeros(len(O.bodies)), label = 'hydroEngine', dead = True),
+	#Measurement, output files
+	PyRunner(command = 'measure()', virtPeriod = 0.1, 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 = 'gravDepo'),
+	#Apply fluid turbulent fluctuations from a discrete random walk model
+	PyRunner(command='turbulentFluctuations()',virtPeriod = 0.01,label = 'turbFluct'),
+	#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='newtonIntegr')
+	]
+#save the initial configuration to be able to recharge the simulation starting configuration easily
+O.saveTmp()
+#run
+#O.run()
+#####################################################################################################################################
+#####################################################  FUNCTION DEFINITION  #########################################################
+#####################################################################################################################################
+
+
+
+######								           ######
+### LET THE TIME FOR THE GRAVITY DEPOSITION AND ACTIVATE THE FLUID AT THE END ###
+######								           ######
+def gravityDeposition(lim):
+	if O.time<lim : return
+	else :		
+		print('\n Gravity deposition finished, apply fluid forces !\n')
+		newtonIntegr.damping = 0.0	# Set the artificial numerical damping to zero
+	   	gravDepo.dead = True	# Remove the present engine for the following
+		hydroEngine.dead = False	# Activate the HydroForceEngine
+		turbFluct.dead = False	# Activate the turbulent fluctuations model
+	return
+###############
+#########################################
+
+
+######							       		           ######
+###   APPLY TURBULENT FLUID VELOCITY FLUCTUATIONS FROM A DISCRETE RANDOM WALK MODEL   ###
+# Associate a random fluctuation to each particles from the value of the Reynolds 	#
+# 	stresses Rxz imposed in simplifiedReynoldsStress = Rxz/densFluid		#
+#Details in the C++ code: function turbulentFluctuationsFluidizedBed in ForceEngine.cpp #
+######								         	   ######
+hydroEngine.simplifiedReynoldStresses = np.ones(ndimz)*4.
+def turbulentFluctuations():
+	hydroEngine.turbulentFluctuationFluidizedBed()
+
+