use of Law2_ScGeom_CapillaryPhys_Capillarity with variable pressure does not change results

Hi,

If I remember correctly, the biggest difference between Law2_ScGeom_CapillaryPhys_Capillarity (which I'm more used to, personnally) and Law2_ScGeom_CapillaryPhys_Capillarity1 is about the underlying capillary files

- Law2_ScGeom_CapillaryPhys_Capillarity1 relies on a unique capillary file, while the other relies on several of them (for different radii)

- the interpolation used to go from the discrete set of data (in the capillary file(s)) to any YADE interaction is also different. The one of Law2_ScGeom_CapillaryPhys_Capillarity1 relies on a Delaunay triangulation, while the other one is a simple (faster ? from my experience) multi-linear one.

Anyway, I would first try much simpler configurations (a liquid bridge between 2 spheres, like [*]) : on that configuration, I'm quite sure you will get different capillary forces for different capillary pressures.

I'm actually expecting your current observations to be quite specific to your current (more complex) configuration.

[*] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/capillaryLaplaceYoung/capillaryBridge.py

Hi Jerôme,

my engine is exactly the same as in the reference [*]:
 #for linear model only
        O.engines=[
                ForceResetter(),
                InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()], avoidSelfInteractionMask = mask2),
                InteractionLoop(
                        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
                        [Ip2_FrictMat_FrictMat_CapillaryPhys()], #for linear model only
                        [Law2_ScGeom_FrictPhys_CundallStrack()], #for linear model only
                ),
                Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=Cap_Pressure),#for linear model only
                NewtonIntegrator(damping=local_damping,gravity=(0,-gravity_y,-gravity_z),label='Newton_integrator'),
                # call the checkUnbalanced function (defined below) every 2 seconds
                PyRunner(command='checkUnbalanced()',iterPeriod=1000,label='checker'),
        ]

I guess that the computeThings() is only for obtaining intermediate output? and does not play a rol ein the fcap calculation?

In my problem, my main concern is the transiont between the packing preparation where the capillary pressure is set to zero in order to have identical packing structure then set to the defined Cap_Pressure.
In my guess, this will not create new contact so no new need of fcap caluation?

maybe I am wrong somewhere?

O.engines does not tell much in itself here. Particle size considerations are at least as much important: with "big" particles, the behavior is anyway insensitive to capillary forces, and I would expect the same result between (dry), (use of Law2_ScGeom_CapillaryPhys_Capillarity with uc1), or (use of Law2_ScGeom_CapillaryPhys_Capillarity with uc2 > uc1)

Regarding your workflow and the sudden activation of capillary forces: usually when you introduce new forces in a discrete system, things can happen (e.g. contacts may create).

As for the "need of fcap calculation", Law2_ScGeom_CapillaryPhys_Capillarity() computes and applies capillary forces at capillary-bonded interactions each time the engine is executed (at each DEM iteration, by default). Not sure I understood what you mean here, actually.

PS : Your reference [*] in #2 can not be the same as mine in #1, it does not help the discussion ;-)

Hi Jérôme,
here the full code and the table for the batch test with only the capillray pressure changes : 1000 and 5000
exactly the same particle positions (611 moving) for all of them after 1913000 steps !!
if I am testing type = 3 it is alaso the same thing : no difference??
of course the results are different between the three cases (1 dry, 2 hertz,3 linear)

So I don't catch why these engine loops are not playing well?

------------------------------------------------
table :
model_type packing_fraction pas_box theta_max converg_min cover_pack_fraction init_seed friction ratio Cap_Pres_Imposed
2 .55 .1 80.0 .002 0.3 625 0.5 2. 1000
2 .55 .1 80.0 .002 0.3 625 0.5 2. 5000

and the code :
# avec directory : +'/'+
# The components of the batch are:
# 1. table with parameters, one set of parameters per line (ccc.table)
# 2. readParamsFromTable which reads respective line from the parameter file
# 3. the simulation muse be run using yade-batch, not yade
#
# $ yade-batch --job-threads=1 03-oedometric-test.table 03-oedometric-test.py
#

# load parameters from file if run in batch
# default values are used if not run from batch
# gravity deposition in box, showing how to plot and save history of data,
# and how to control the simulation while it is running by calling
# python functions from within the simulation loop

#pas_box theta_max converg_min cover_pack_fraction init_seed friction ratio
#model_type bottom_cover pas_box theta_max converg_min cover_pack_fraction init_seed friction ratio
#3 .55 .1 80.0 .002 0.3 625 0.5 2.5 5000
readParamsFromTable(model_type=1, packing_fraction=0.6,pas_box=0.1, theta_max=45.0, converg_min=0.00005,cover_pack_fraction=.2,init_seed=10, friction=0.5,ratio=2.5,Cap_Pres_Imposed=1000)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *

# import yade modules that we will use below
from yade import pack, plot, export,math
global ratio,nombre,friction # size ratio between the glued spheres and the moving ones
global cover_pack_fraction ,converg_min, init_seed,nombre_moving # coverage percent for moving spheres
global i_pas, box_size,pas_box,rayon,step0,step1,step2,step_precedent,gravity_y,gravity_z,theta_max,theta_max,nom_file,filename_yade,traitement_file,str_Angle
global angle_actif,angle

#some parameters passed by batch_table:

####batch : friction = 0.5
####batch : theta_max = 30.0
####batch : pas_box = 0.1
####batch : ratio = 3 # size ratio between the glued spheres and the moving ones
####batch : cover_pack_fraction = 0.2 # coverage percent for moving spheres
####batch : init_seed=10
####batch : packing_fraction = 70./100.
####batch : model_type = 1
#some parameters:
shear_modulus = 1e5
poisson_ratio = 0.3
young_modulus = 2*shear_modulus*(1+poisson_ratio)
local_damping = 0.01
viscous_normal = 0.021
viscous_shear = 0.8*viscous_normal
position_init = numpy.arange(1200*4, dtype=float)
position_init = position_init.reshape((1200,4))
numpy.zeros_like(position_init)
Cap_Pressure = 0.0 #valeur initiale
angle = math.atan(friction) # value in radians for friction Angle
angle_init = math.atan(0.10)
angle_actif = angle_init
# initialisation coordonnees initiales
#newTable("position_init",1200,4) # Create a new table with 5 rows and 3 column
#creating a material (FrictMat):
#frictionAngle(=.5)
# Contact friction angle (in radians). Hint : use ‘radians(degreesValue)’ in python scripts.

id_SphereMat=O.materials.append(FrictMat(young=young_modulus,poisson=poisson_ratio,density=2500,frictionAngle=angle_actif,label="glass_beads"))
SphereMat=O.materials[id_SphereMat]

#1 = >hertz model sans capillary forces
if model_type == 1:#hertz model with capillary forces
 converg_min = converg_min
elif model_type == 2: #hertz model with capillary forces
 converg_min = 10.0*converg_min
else: #for linear model only
 converg_min = 10*converg_min
box_size = 4.0
i_pas = 0

step0 = 0
step1 = 0
step2 = 0
step_precedent = 0
mask1=0b01
mask2=0b10
mask3=0b11

gravity_y = 9.81*sin( pas_box*i_pas*3.14/180.0)
gravity_z = 9.81*cos( pas_box*i_pas*3.14/180.0)

O.reset

# create rectangular box from facets
O.bodies.append(geom.facetBox((box_size/2.0,box_size/2.0,box_size*0.62),(box_size/2.0,box_size/2.0,box_size*0.62),wallMask=31))

# create empty sphere packing
# sphere packing is not equivalent to particles in simulation, it contains only the pure geometry

rayon = 0.025
height = 5. * rayon
nombre= int(packing_fraction*(box_size*box_size)/(rayon*rayon*3.14159) )
print("nombre = ",nombre)
print("friction angle = ",angle_actif)
print("Model type (1 sec- 2 hertz, 3 lineaire) = ",model_type)
sp=pack.SpherePack()
# generate randomly spheres with uniform radius distribution
sp.makeCloud((0,0,0),(box_size,box_size,height),rMean=rayon,rRelFuzz=.0005,num= nombre,seed=init_seed)

# add the sphere pack to the simulation
sp.toSimulation(color=(1,0,0),mask=mask1,material=SphereMat)

# simulation loop
#define engines:
if model_type == 1:
 O.engines=[
   ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()], avoidSelfInteractionMask = mask2),
  InteractionLoop(
   # handle sphere+sphere and facet+sphere collisions
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_FrictPhys(frictAngle=angle_actif)],
                        [Law2_ScGeom_FrictPhys_CundallStrack()]
# [Law2_L3Geom_FrictPhys_ElPerfPl()]
  ),
  NewtonIntegrator(damping=local_damping ,gravity=(0,-gravity_y,-gravity_z),label='Newton_integrator'),
  # call the checkUnbalanced function (defined below) every 2 seconds
  PyRunner(command='checkUnbalanced()',iterPeriod=1000,label='checker'),
 ]
elif model_type == 2: #hertz model with capillary forces
 O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()], avoidSelfInteractionMask = mask2),
  InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_MindlinCapillaryPhys(label='ContactModel')],#for hertz model only
   [Law2_ScGeom_MindlinPhys_Mindlin()] #for hertz model only
  ),
  Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=Cap_Pressure ),
                #It is also necessary to disable interactions removal by the constitutive law (Law2).
                #The only combinations of laws supported are currently capillary law + Law2_ScGeom_FrictPhys_CundallStrack
                #and capillary law + Law2_ScGeom_MindlinPhys_Mindlin (and the other variants of Hertz-Mindlin).
  #for hertz model only Value of the capillary pressure Uc defined as Uc=Ugas-Uliquid 10000 initial
  #Capillary pressure uc, to be defined equal to Law2_ScGeom_CapillaryPhys_Capillarity.capillaryPressure.
  NewtonIntegrator(damping=local_damping ,gravity=(0,-gravity_y,-gravity_z),label='Newton_integrator'),
  # call the checkUnbalanced function (defined below) every 2 seconds
  PyRunner(command='checkUnbalanced()',iterPeriod=1000,label='checker'),
 ]
 ContactModel.betan=viscous_normal
 ContactModel.betas=viscous_shear
 ContactModel.useDamping=True
else: #for linear model only
 O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()], avoidSelfInteractionMask = mask2),
  InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_CapillaryPhys()], #for linear model only
   [Law2_ScGeom_FrictPhys_CundallStrack()], #for linear model only
  ),
  Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=Cap_Pressure),#for linear model only
  NewtonIntegrator(damping=local_damping,gravity=(0,-gravity_y,-gravity_z),label='Newton_integrator'),
     # call the checkUnbalanced function (defined below) every 2 seconds
  PyRunner(command='checkUnbalanced()',iterPeriod=1000,label='checker'),
 ]

# timesteps
O.dt=.5*PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy=True

nom_file=str(model_type)+"_"+str(packing_fraction)+"_"+str(cover_pack_fraction)+"_"+str(ratio)+'_'+str(friction)+'_'+str(pas_box)+'_'+str(converg_min)+'_'+str(init_seed)+'_'+str(Cap_Pres_Imposed)
filename_yade=nom_file+'.yade'
traitement_file = nom_file+'.traitement'
traitement = open(traitement_file, 'w')
traitement.close()
# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
 global converg_min,friction,theta_max,step0,step1,step2, pas_box,i_pas,step_precedent,rayon,init_seed,ratio,cover_pack_fraction,nom_file,filename_yade,nombre_moving
# print("iteration= ",utils.unbalancedForce(),O.iter,step_precedent)# numéro de l'étape
 if (O.iter)<(5000): return
 if (O.iter-step0)<(5000): return
 if (O.iter-step_precedent)>100000 and step_precedent>0 and step2==1:
  O.save('3rd-step_interrompu_'+filename_yade)
  O.pause()
#########################################################################
# the rest will be run only if unbalanced is < <converg_min (stabilized packing)
 if utils.unbalancedForce()>converg_min: return
# print ("bal iter steps :",utils.unbalancedForce(),O.iter,step0,step1,step2)
#########################################################################
 if step0==0:
  step0 = O.iter
  print ("len-O.bodies = ",len(O.bodies))
  for b in O.bodies:
   #b1.state.vel = Vector3(0,0,0)
   b.state.vel[0] = 0 # stop
   b.state.vel[1] = 0 # tous les
   b.state.vel[2] = 0 # deplacements translations
   b.state.angVel = Vector3(0,0,0)# stop les rotations
   #b.angVel[0]=0
   #b.angVel[1]=0
   #b.angVel[2]=0
   b.state.blockedDOFs='xyzXYZ'# réinisialise les accelerations
   b.state.dynamic=False #réinisialise les vitesses
  O.save('1st-step_'+filename_yade)
# save file text mode; beginning of the run
# export.textExt('./output_fin.txt',format='x_y_z_r_attrs',attrs=('b.state.pos.norm()','b.state.pos'),comment='dstN dstV_x dstV_y dstV_z')
# export.textExt('./output_fin.txt',format='x_y_z_r_attrs',attrs=['b.state.vel[0]','b.state.vel[1]','b.state.vel[2]'],comment='Vx Vy Vz')
  print ("fin de step 1: ",step0)
##################creation deuxieme couche #######################################################
 if step0 ==O.iter and step1==0:
  step1=1
  nombre_moving = int(cover_pack_fraction*(box_size*box_size)/(ratio*rayon*ratio*rayon*3.14159))
  print ("nombre_moving =", nombre_moving)
#generate particles:
  sp=pack.SpherePack()
# generate randomly spheres with uniform radius distribution
  sp.makeCloud((rayon,rayon,2.*rayon),(box_size-rayon,box_size-rayon,(3.2+.1)*rayon*ratio),rMean=rayon*ratio,rRelFuzz=.005*ratio,num= nombre_moving,seed=2*init_seed)
# O.bodies.append([sphere(c,r,material=SphereMat,color=(0,2,0)) for c,r in sp])
# add the sphere pack to the simulation
  sp.toSimulation(color=(0,1,0),mask = mask3,material=SphereMat)#,material=SphereMat
# export.textExt('./output_creation_2packs.txt',format='x_y_z_r_attrs',attrs=['b.state.vel[0]','b.state.vel[1]','b.state.vel[2]'],comment='V_x V_y V_z')
  print ("len-O.bodies (2 packs) = ",len(O.bodies))

#########################################################################
 if O.iter>step0 and step1==1:
   step1=O.iter
   step2=1
   step_precedent=O.iter
   O.save('2nd-step_'+filename_yade)
   angle_actif = angle
   O.materials[0].frictionAngle=angle_actif # radians
   # for existing contacts, set contact friction directly
   for i in O.interactions: i.phys.tangensOfFrictionAngle=tan(angle_actif)
   print("friction angle (step2) = ",angle_actif)
   pos_init()
# save file text mode; beginning of the run
# export.textExt('./output_fin.txt',format='x_y_z_r_attrs',attrs=('b.state.pos.norm()','b.state.pos'),comment='dstN dstV_x dstV_y dstV_z')
# export.textExt('./output_fin-2ndstep.txt',format='x_y_z_r_attrs',attrs=['b.state.vel[0]','b.state.vel[1]','b.state.vel[2]'],comment='Vx Vy Vz')
   print ("fin de step2 : ", step1)
   Cap_Pressure=Cap_Pres_Imposed
#########################################################################
 if O.iter>step1 and step2==1:
  if O.iter<step_precedent+10000 : return
  step_precedent=O.iter
  num_Angle=pas_box*i_pas
  str_Angle= f"{(num_Angle):.3f}"
  pos_moved()
# filename1='./output_'+nom_file+'_'+str(pas_box*i_pas)+'_'+str(i_pas)+'.txt'
  filename1='./output_'+nom_file+'_'+str_Angle+'_'+str(i_pas)+'.txt'
  filenamevtk='./output_'+nom_file+'_'+str(i_pas)+'.vtk'
  export.textExt(filename1,format='x_y_z_r_attrs',attrs=['b.state.vel[0]','b.state.vel[1]','b.state.vel[2]'],comment='V_x V_y V_z')
# export.text('/tmp/test.txt')
# # now open paraview, click "Tools" menu -> "Python Shell", click "Run Script", choose "pv_section.py" from this directiory
# or just run python pv_section.py
# see /home/oger/yade_full_install_local/trunk/examples/test/paraview-spheres-solid-section/pv_section.py
# text2vtk and text2vtkSection function can be copy-pasted from yade/py/export.py into separate py file to avoid executing yade or to use pure python
  export.text2vtk(filename1,filenamevtk)
  i_pas = i_pas+ 1
  gravity_y = 9.81*sin( pas_box*i_pas*3.14/180.0)
  gravity_z = 9.81*cos( pas_box*i_pas*3.14/180.0)
  Newton_integrator.gravity=Vector3(0,-gravity_y,-gravity_z)
  print ("info : ", (pas_box*i_pas), gravity_y, gravity_z, O.iter)
  ForceResetter()
  #print("iteration= ", O.iter)# numéro de l'étape
  if (pas_box*i_pas)>=theta_max :
   # save output file txt after end of second sedimentation
####batch : export.textExt('./output_fin_3rd-step.txt',format='x_y_z_r_attrs',attrs=['b.state.vel[0]','b.state.vel[1]','b.state.vel[2]'],comment='V_x V_y V_z')
   O.save('3rd-step_'+filename_yade)
   O.pause()

def pos_init():
 global position_init,nombre_mobile,nombre_moving
 nombre_mobile=0
 for b in O.bodies:
  if b.mask==mask3:
   position_init[nombre_mobile][0]=b.state.pos[0]
   position_init[nombre_mobile][1]=b.state.pos[1]
   position_init[nombre_mobile][2]=b.state.pos[2]
   position_init[nombre_mobile][3]=b.shape.radius
# print( position_init[nombre_mobile][0], position_init[nombre_mobile][1], position_init[nombre_mobile][2], position_init[nombre_mobile][3])
   nombre_mobile = nombre_mobile+1

def pos_moved():
 global position_init,str_Angle,nombre_moving,traitement_file,i_pas,ratio,pas_box
 nombre_mobile = 0
 nombre_moved = 0
 for b in O.bodies:
  if b.mask==mask3:
   dist=((position_init[nombre_mobile][0]-b.state.pos[0])*(position_init[nombre_mobile][0]-b.state.pos[0]) + \
   (position_init[nombre_mobile][1]-b.state.pos[1])*(position_init[nombre_mobile][1]-b.state.pos[1]) + \
   (position_init[nombre_mobile][2]-b.state.pos[2])*(position_init[nombre_mobile][2]-b.state.pos[2]))
   if dist>( position_init[nombre_mobile][3]* position_init[nombre_mobile][3]):
    nombre_moved=nombre_moved + 1
# print("====>",nombre_mobile,nombre_moved)
   nombre_mobile = nombre_mobile+1
 traitement = open(traitement_file, 'a')
 print (i_pas,nombre_moved,nombre_moving,i_pas*pas_box,nombre_moved/nombre_moving, file=traitement)
 traitement.close()
# save tmp pour rerun convergence
O.saveTmp()

O.run()
# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

i still need an answer so:
please help me by just checking if the series of mechanical parameters attribution are correct as they are not playing the exact rules right now!

This question was expired because it remained in the 'Open' state without activity for the last 15 days.