yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #12844
[Branch ~yade-pkg/yade/git-trunk] Rev 3935: Capillary scripts commit (http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg12103.html)...
------------------------------------------------------------
revno: 3935
committer: jduriez <jerome.duriez@xxxxxxxxxxx>
timestamp: Wed 2016-09-21 10:48:23 -0600
message:
Capillary scripts commit (http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg12103.html). Located in a new examples/capillaryLaplaceYoung folder intended to illustrates capillary Law2. With one new simulation script example and another old one, moved here
removed:
examples/CapillaryPhys-example.py
added:
examples/capillaryLaplaceYoung/
examples/capillaryLaplaceYoung/CapillaryPhys-example.py
examples/capillaryLaplaceYoung/README.txt
examples/capillaryLaplaceYoung/capillaryBridge.py
examples/capillaryLaplaceYoung/solveLaplace_uc.m
examples/capillaryLaplaceYoung/solveLiqBridge.m
examples/capillaryLaplaceYoung/writesCapFile.m
--
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
=== removed file 'examples/CapillaryPhys-example.py'
--- examples/CapillaryPhys-example.py 2013-05-30 20:17:45 +0000
+++ examples/CapillaryPhys-example.py 1970-01-01 00:00:00 +0000
@@ -1,74 +0,0 @@
-#!/usr/bin/python
-# -*- coding: utf-8 -*-
-
-"""
-To run this script you need to have all 10 text files from https://yade-dem.org/wiki/CapillaryTriaxialTest
-in the same folder as you run this script in console!
-
-This script shows how to use Law2_ScGeom_CapillaryPhys_Capillarity. The user can switch between hertz and
-linear model by setting "model_type"
-"""
-
-model_type = 1 #1=Hertz model with capillary forces, else linear model with capillary model
-
-#some parameters:
-shear_modulus = 1e5
-poisson_ratio = 0.3
-young_modulus = 2*shear_modulus*(1+poisson_ratio)
-friction = 0.5
-angle = atan(friction)
-local_damping = 0.01
-viscous_normal = 0.021
-viscous_shear = 0.8*viscous_normal
-lowercorner = Vector3(0,0,0)
-uppercorner = Vector3(0.002,0.002,0.004)
-
-#creating a material (FrictMat):
-id_SphereMat=O.materials.append(FrictMat(young=young_modulus,poisson=poisson_ratio,density=2500,frictionAngle=angle))
-SphereMat=O.materials[id_SphereMat]
-
-#generate particles:
-sp=pack.SpherePack()
-sp.makeCloud(lowercorner,uppercorner,.0002,rRelFuzz=.3)
-O.bodies.append([sphere(c,r,material=SphereMat) for c,r in sp])
-
-#generate boundary:
-O.bodies.append(geom.facetBox(uppercorner/2,uppercorner/2,wire=True,fixed=True,material=SphereMat))
-
-#define engines:
-if model_type == 1:#hertz model with capillary forces
- O.engines=[
- ForceResetter(),
- InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
- 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=10000),#for hertz model only
- NewtonIntegrator(damping=local_damping,gravity=(0,0,-9.81)),
- ]
- ContactModel.betan=viscous_normal
- ContactModel.betas=viscous_shear
- ContactModel.useDamping=True
-else:
- O.engines=[
- ForceResetter(),
- InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
- 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=10000),#for linear model only
- NewtonIntegrator(damping=local_damping,gravity=(0,0,-9.81)),
- ]
-
-#set time step and run simulation:
-O.dt=0.5*PWaveTimeStep()
-
-from yade import qt
-qt.View()
-print('Press PLAY button')
-#O.run(10000,True)
-
=== added directory 'examples/capillaryLaplaceYoung'
=== added file 'examples/capillaryLaplaceYoung/CapillaryPhys-example.py'
--- examples/capillaryLaplaceYoung/CapillaryPhys-example.py 1970-01-01 00:00:00 +0000
+++ examples/capillaryLaplaceYoung/CapillaryPhys-example.py 2016-09-21 16:48:23 +0000
@@ -0,0 +1,74 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+"""
+To run this script you need to have all 10 text files from https://yade-dem.org/wiki/CapillaryTriaxialTest
+in the same folder as you run this script in console!
+
+This script shows how to use Law2_ScGeom_CapillaryPhys_Capillarity. The user can switch between hertz and
+linear model by setting "model_type"
+"""
+
+model_type = 1 #1=Hertz model with capillary forces, else linear model with capillary model
+
+#some parameters:
+shear_modulus = 1e5
+poisson_ratio = 0.3
+young_modulus = 2*shear_modulus*(1+poisson_ratio)
+friction = 0.5
+angle = atan(friction)
+local_damping = 0.01
+viscous_normal = 0.021
+viscous_shear = 0.8*viscous_normal
+lowercorner = Vector3(0,0,0)
+uppercorner = Vector3(0.002,0.002,0.004)
+
+#creating a material (FrictMat):
+id_SphereMat=O.materials.append(FrictMat(young=young_modulus,poisson=poisson_ratio,density=2500,frictionAngle=angle))
+SphereMat=O.materials[id_SphereMat]
+
+#generate particles:
+sp=pack.SpherePack()
+sp.makeCloud(lowercorner,uppercorner,.0002,rRelFuzz=.3)
+O.bodies.append([sphere(c,r,material=SphereMat) for c,r in sp])
+
+#generate boundary:
+O.bodies.append(geom.facetBox(uppercorner/2,uppercorner/2,wire=True,fixed=True,material=SphereMat))
+
+#define engines:
+if model_type == 1:#hertz model with capillary forces
+ O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
+ 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=10000),#for hertz model only
+ NewtonIntegrator(damping=local_damping,gravity=(0,0,-9.81)),
+ ]
+ ContactModel.betan=viscous_normal
+ ContactModel.betas=viscous_shear
+ ContactModel.useDamping=True
+else:
+ O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
+ 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=10000),#for linear model only
+ NewtonIntegrator(damping=local_damping,gravity=(0,0,-9.81)),
+ ]
+
+#set time step and run simulation:
+O.dt=0.5*PWaveTimeStep()
+
+from yade import qt
+qt.View()
+print('Press PLAY button')
+#O.run(10000,True)
+
=== added file 'examples/capillaryLaplaceYoung/README.txt'
--- examples/capillaryLaplaceYoung/README.txt 1970-01-01 00:00:00 +0000
+++ examples/capillaryLaplaceYoung/README.txt 2016-09-21 16:48:23 +0000
@@ -0,0 +1,35 @@
+This folder illustrates the use of Law2_ScGeom_CapillaryPhys_Capillarity engine to describe capillary interaction between particle pairs in YADE (see also the wiki page https://www.yade-dem.org/wiki/CapillaryTriaxialTest)
+
+
+I. Simulation scripts examples
+------------------------------
+
+Two examples of simulations using Law2_ScGeom_CapillaryPhys_Capillarity are herein provided:
+
+- CapillaryPhys-example.py illustrates the mutual attraction due to capillary interaction in a simple packing
+- capillaryBridge.py defines and let evolve a capillary bridge between two particles only
+
+
+II. Capillary files scripts
+---------------------------
+
+As explained on the wiki page (see above) and in Law2_ScGeom_CapillaryPhys_Capillarity doc, so-called capillary files are required in order to use Law2_ScGeom_CapillaryPhys_Capillarity engine. These capillary files include capillary bridges configurations data (for positiv capillary pressure values). A version of capillary files can be downloaded from the wiki page, they consider a zero contact (wetting) angle and some given particle radius ratio.
+
+Also, three .m files are herein provided in order to enable users to build their own capillary files for any contact angle or particle radius ratio, or to study capillary bridges such as such (solveLiqBridge.m). A capillary files generation requires launching writesCapFile() (in writesCapFile.m) only, after providing a couple of parameters therein (see description of writesCapFile.m below). The two other files (solveLiqBridge.m and solveLaplace_uc.m) define functions used by writesCapFile() and do not require any user input. Specifically, the roles of the three .m files are as follows:
+
+
+- solveLiqBridge.m solves the Laplace-Young equation for one given bridge, defined in terms of the input attributes of the solveLiqBridge function (see therein). The solveLiqBridge function is usually called by other files (see below) during capillary files generation, however it can also be executed on its own in order to study (e.g. plot) capillary bridge profile.
+
+Code comments include references to:
+ * Duriez2016: J. Duriez and R. Wan, Contact angle mechanical influence for wet granular soils, currently under Review in Acta Geotechnica
+ * Lian1993: G. Lian and C. Thornton and M. J. Adams, A Theoretical Study of the Liquid Bridge Forces between Two Rigid Spherical Bodies, Journal of Colloid and Interface Science, 161(1), 1993
+ * Scholtes2008 (french): L. Scholtes, Modelisation Micro-Mecanique des Milieux Granulaires Partiellement Satures, PhD Thesis from Institut polytechnique de Grenoble, 2008
+ * Soulie2005 (french): F. Soulie, Cohesion par capillarite et comportement mecanique de milieux granulaires, PhD Thesis from Universite Montpellier II, 2005
+ * Soulie2006: F. Soulie and F. Cherblanc and M. S. El Youssoufi and C. Saix, Influence of liquid bridges on the mechanical behaviour of polydisperse granular materials, International Journal for Numerical and Analytical Methods in Geomechanics, 30(3), 2006
+
+
+- solveLaplace_uc.m calls several times solveLiqBridge in order to compute all possible bridges configurations for given contact angle, capillary pressure, and particle radius ratio. It is usually called by writesCapFile (see below) during capillary files generation, however it can also be executed on its own. In particular, it may output text files including bridges data for one given capillary pressure.
+
+
+- writesCapFile.m, finally, is the conductor during the capillary files generation, and is the only one that actually requires the user's attention for such task. Parameters are to be defined by the user directly in the .m file (from l. 10 to 30) before launching for data files generation.
+ Generating a complete set of 10 capillary files typically requires few days in terms of computation time with MATLAB and much (~ 7*) more with Octave. The .m files makes wide use of iterative procedure and a greater performance for MATLAB is indeed expected. In case Octave is still chosen for generation, you may want uncomment some Octave specific command ('fflush(..)') in writesCapFile.m in order for the terminal to update you on the progress.
=== added file 'examples/capillaryLaplaceYoung/capillaryBridge.py'
--- examples/capillaryLaplaceYoung/capillaryBridge.py 1970-01-01 00:00:00 +0000
+++ examples/capillaryLaplaceYoung/capillaryBridge.py 2016-09-21 16:48:23 +0000
@@ -0,0 +1,50 @@
+# capillary bridge example between 2 spheres
+
+
+r1,r2 = 1e-4,4e-4 # better take small particles
+z1,z2=0,0.95*(r1+r2)
+
+O.bodies.append(sphere(center=Vector3(0,0,z1),radius=r1,dynamic=0))
+O.bodies.append(sphere(center=Vector3(0,0,z2),radius=r2,dynamic=0))
+
+O.engines=[ForceResetter()
+ ,InsertionSortCollider([Bo1_Sphere_Aabb()])
+ ,InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom()],
+ [Ip2_FrictMat_FrictMat_CapillaryPhys()],
+ [Law2_ScGeom_FrictPhys_CundallStrack(neverErase=1)]
+ )
+ ,Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=1e3)
+ ,NewtonIntegrator()
+ ,GlobalStiffnessTimeStepper()
+ ,PyRunner(command='computeThings()',iterPeriod=1)
+ ]
+
+from yade import plot
+def computeThings():
+ if O.interactions.has(0,1) and O.interactions[0,1].isReal:
+ i = O.interactions[0,1]
+ un = i.geom.penetrationDepth
+ vM = i.phys.vMeniscus
+ delta1,delta2 = i.phys.Delta1,i.phys.Delta2
+ fCap = i.phys.fCap.norm()
+ nn11,nn33 = i.phys.nn11,i.phys.nn33
+ else:
+ un,vM,delta1,delta2,fCap,nn11,nn33 = float('nan'),float('nan'),float('nan'),float('nan'),float('nan'),float('nan'),float('nan')
+
+ plot.addData(delta1 = delta1, delta2 = delta2, fCap = fCap
+ , it =O.iter, un = un ,vM = vM
+ , nn11 = nn11, nn33 = nn33)
+
+plot.plots={'un':'fCap','un ':'vM',' un':('delta1','delta2'),' un ':('nn11','nn33')}
+plot.plot()
+
+# spheres get closer:
+O.bodies[1].state.vel=Vector3(0,0,-0.002)
+O.run(500,wait=1)
+
+# spheres move apart:
+O.bodies[1].state.vel=Vector3(0,0,0.02)
+O.run(4000,wait=1)
+
+
=== added file 'examples/capillaryLaplaceYoung/solveLaplace_uc.m'
--- examples/capillaryLaplaceYoung/solveLaplace_uc.m 1970-01-01 00:00:00 +0000
+++ examples/capillaryLaplaceYoung/solveLaplace_uc.m 2016-09-21 16:48:23 +0000
@@ -0,0 +1,105 @@
+% J. Duriez (jerome.duriez@xxxxxxxxxxx)
+
+function stabSol = solveLaplace_uc(theta,rRatio,uStar,delta1Cons,deltaZ,save)
+% Solves Laplace-Young equation for given r, uc* and theta
+% Making recursively use of solveLiqBridge function, the present function
+% identifies all possible bridges configurations for given r, uc*, and
+% theta
+
+% Input attributes (see also solveLiqBridge function):
+% theta: contact angle (deg)
+% rRatio: rMax / rMin for the considered particles pair
+% uStar: dimensionless capillary pressure
+% delta1Cons: vector of delta1 (see solveLiqBridge) values to consider
+% deltaZ: see solveLiqBridge
+% save: 0/1 to save text files including the capillary bridges data (if 1)
+
+% Returns a table (matrix) of capillary bridges data (see sol_u below),
+% possibly (in case save=1) also written in text files.
+
+
+
+% Removal of possible d1=0 since such bridge is physically impossible, and
+% d1 = 0 might make trouble if theta = 0 (leads to inf. initial rhoPrime):
+delta1Cons = delta1Cons(delta1Cons~=0);
+% And of d1 values such that d1 + theta > 90
+delta1Cons = delta1Cons(delta1Cons + theta <= 90);
+
+
+%------------ All possible bridges for these r, uc*, theta: ---------------
+
+% data table initialization:
+sol_u = -ones(length(delta1Cons),9);
+% 9 col. are [dist uc(=cst) V F delta1 delta2 E nn11 nn33]
+
+
+minDzUsed = deltaZ;
+% tic
+for i=1:length(delta1Cons)
+ dzUsed = deltaZ;
+ delta1 = delta1Cons(i);
+ [dist,vol,force,delta1,delta2,eStar,nn11,nn33,out] = solveLiqBridge(rRatio,theta,uStar,delta1,dzUsed,0,0,0);
+ while out == 0
+ dzUsed = dzUsed / 2;minDzUsed = min(dzUsed,minDzUsed);
+% disp(['Refining deltaZ to ',num2str(dzUsed),' in solveLaplace'])
+ [dist,vol,force,delta1,delta2,eStar,nn11,nn33,out] = solveLiqBridge(rRatio,theta,uStar,delta1,dzUsed,0,0,0);
+ end
+ sol_u(i,:) = [dist,uStar,vol,force,delta1,delta2,eStar,nn11,nn33];
+
+end
+% toc
+
+% Get rid of unphysical solutions with negativ distances and/or volume:
+physSol = sol_u( sol_u(:,1)>=0 & sol_u(:,3)>=0,: ); % NB: vol>=0 is true for vol = Inf
+
+% Get rid of diverged solutions with infinite volume (should <=> infinite
+% profile):
+nonDvSol = physSol(isfinite(physSol(:,3)),:);
+% nNonDvSol = length(nonDvSol(:,1));
+% disp(['We had to suppress ',num2str((1-nNonDvSol/nPhysSol)*100),' % of diverged "solutions"'])
+
+% Get rid of unstable physical solutions:(those with biggest volumes)
+% We use volume values for this purpose, see e.g. Duriez2016
+distRupt = max(nonDvSol(:,1));
+% eRupt = nonDvSol(nonDvSol(:,1)==distRupt,7); % does not work since eRupt can be a global maximum
+% (when two branches are increasing with d*)
+vRupt = nonDvSol(nonDvSol(:,1)==distRupt,3);
+stabSol = nonDvSol(nonDvSol(:,3)>=vRupt,:);
+
+
+% ----------- Save of data in text files -----------
+if save ~=0
+% three text files are generated: one, "Raw", includes all computed data;
+% another, "Phys", restricts to physical solutions (with positiv distance
+% and volume); the last one, "Stab" restricts to the stable physical
+% solutions
+ if minDzUsed == deltaZ
+ strDz = ['With deltaZ = ',num2str(deltaZ)];
+ else
+ strDz = ['With deltaZ between ',num2str(deltaZ),' and ',num2str(minDzUsed)];
+ end
+ nomFic = ['capDataRaw_r',num2str(rRatio),'_theta',num2str(theta),'_ucStar',num2str(uStar),'.txt'];
+ head_wE = ['dist*\tuc*\tV*\tF*\tdelta1 (deg)\tdelta2(deg)\tE* (-)\tn11 (-)\tn33 (-)\t',strDz,'\n'];
+ writeFile(nomFic,head_wE,sol_u)
+
+ nomFic = ['capDataPhys_r',num2str(rRatio),'_theta',num2str(theta),'_ucStar',num2str(uStar),'.txt'];
+ writeFile(nomFic,head_wE,nonDvSol)
+
+ nomFic = ['capDataStab_r',num2str(rRatio),'_theta',num2str(theta),'_ucStar',num2str(uStar),'.txt'];
+ head_woE = ['dist*\tuc*\tV*\tF*\tdelta1 (deg)\tdelta2(deg)\tn11 (-)\tn33 (-)\t',strDz,'\n'];
+ writeFile(nomFic,head_woE,stabSol(:,[1:6,8:9]))
+end
+% ---------------------------------------------------
+
+end
+
+
+function writeFile(nomFic,headLine,data)
+
+fic = fopen(nomFic,'w');
+fprintf(fic,headLine);
+fclose(fic);
+dlmwrite(nomFic,data,'-append','delimiter','\t')
+
+end
+
=== added file 'examples/capillaryLaplaceYoung/solveLiqBridge.m'
--- examples/capillaryLaplaceYoung/solveLiqBridge.m 1970-01-01 00:00:00 +0000
+++ examples/capillaryLaplaceYoung/solveLiqBridge.m 2016-09-21 16:48:23 +0000
@@ -0,0 +1,279 @@
+% J. Duriez (jerome.duriez@xxxxxxxxxxx)
+
+function [dist,vol,force,delta1,delta2,eStar,nn11,nn33,out] = solveLiqBridge(rRatio,theta,uStar,delta1,deltaZ,plot,speak,plot3D)
+%Computes one single liquid (capillary) bridge
+% Lian's method: iterative determination using Taylor expansion with
+% finite difference (FD) increment = deltaZ (a function attribute)
+
+% Attributes: rRatio = R2 / R1 >=1, theta = contact angle (deg), uStar =
+% dimensionless cap. pressure (uc*R2/gamma), delta1 = filling angle on
+% smallest particle (deg), deltaZ = increment of FD scheme (e.g. 2e-6)
+% Last ones = booleans to resp. plot the profile in 2D, output messages,
+% and plot the liq bridge in 3D
+
+% Outputs: dist = interparticle distance ;
+% vol = dimensionless volume (= vol% /R2^3)
+% force = dimensionless force = force/(2*pi*gamma*R2)
+% delta1/2 = filling angles on the 2 particles (degrees)
+% eStar = free energy/(gamma*% R2^2)
+% nn11 = 11-term of integral(ni nj dS) / R2^2 = 22-term
+% nn33 = 33-term, with 3 axis = the meniscus orientation
+% NB: surface / R2^2 = 2*nn11 + nn33
+% out = 1 or 0 depending whether things went fine (see below): 1 if yes
+
+% close all
+global cstC
+% cstC is the dimensionless capillary force: constant all along the profile
+% see [9] Lian1993, (6) Duriez2016, or (2.51) Soulie2005 ~ (10) Soulie2006..
+cstC = 1/rRatio * sin(radians(delta1)) * sin(radians(delta1+theta)) + 1/2 * uStar * rRatio^-2 * sin(radians(delta1))^2;
+
+% Use of cstC to get the right filling angle delta2:
+fun = @(delta2)cstC_delta2(delta2,theta,uStar);
+delta2 = fzero(fun,[0 90]); % fzero does not catch multiple roots, if it occurs
+
+% Discrete set of profile values:
+rho = -ones(30000,1); % preallocating 26000 or 30000 makes 90 profile computations take ~ 10.8s instead of ~ 14.5
+% Discrete set of slope values:
+rhoPrime = -ones(30000,1);
+
+step = 1;
+% Boundary conditions on left contact line: see left (smallest) particle
+rho(1) = 1/rRatio * sin(radians(delta1));
+rhoPrime(1) = - 1 / tan(radians(delta1+theta));
+% Boundary condition on right contact line: see right (biggest) particle
+rhoRight = sin(radians(delta2));
+
+%-------------------------------------------------------------------------
+% Finite diff. scheme to compute whole profile ie compute rho and rhoPrime
+% See Lian1993, Duriez2016, etc..
+
+% Some remarks about next loop:
+% - this loop may extend the size of rho without problem
+% - just keep a simple and efficient test to stop computations. A more restrictive
+% one may lead to go beyond the right contact line and lead to divergence..
+while rho(step) < 1.00001*rhoRight
+ rho2d = rhoSecond(rho(step),rhoPrime(step),uStar); % all terms are at "i"
+ rho(step+1) = rho(step) + deltaZ * rhoPrime(step) + 1/2 * deltaZ^2 * rho2d;
+
+ rhoPrime(step+1) = drho(rho(step+1),rho(step),deltaZ,rho2d); % i+1 value of rhoPrime, function of rho_{i+1}, rho_i, and rho2d_i. Used next time in loop.
+ step=step+1;
+end
+%-------------------------------------------------------------------------
+
+% Two following lines to suppress extra data due to preallocation. Keep
+% this order !
+rhoPrime = rhoPrime(rho>0);
+rho = rho(rho>0);
+
+% --------------- Output computations from now on: ------------------------
+
+d1 = radians(delta1); d2 = radians(delta2);
+
+% Inter particle dimensionless distance:
+% see (7) Duriez2016, (33) Scholtes2008, etc. (not (5) Soulie2006..)
+lastZ = deltaZ*(length(rho) -1);
+dist = lastZ - 1/rRatio * ( 1-cos(d1) ) - ( 1-cos(d2) );
+
+% Dimensionless capillary force:
+force = cstC;
+
+% Dimensionless volume:
+% Cf (4) Soulie2006, (34) Appendix Scholtes2008, (8) Duriez2016 etc. :
+vol = pi*sum(rho(2:length(rho)).^2)*deltaZ;
+vol = vol - pi/3 * rRatio^(-3) * ( 1-cos(d1) )^2 * ( 2+cos(d1) );
+vol = vol - pi/3 * ( 1-cos(d2) )^2 * ( 2+cos(d2) );
+
+% Capillary bridge dimensionless free energy:
+% (12) Duriez2016 rather than [35] Lian1993 (was for cst volume stability)
+dArea = rho .* ( 1+rhoPrime.^2).^(1/2); % ~ infinitesimal liquid gas area
+eStar = 2*pi * sum(dArea(2:length(dArea))) * deltaZ + uStar * vol; %+/- uStar changes the values but not the shape
+eStar = eStar - 2*pi*cos(radians(theta)) * ( (1-cos(d1))/rRatio^2 + 1 - cos(d2) );
+
+rhoRed = rho(2:length(rho)); rhoPrRed = rhoPrime(2:length(rho));
+
+% Surface:
+surf = 2*pi * sum( rhoRed .* (1 + rhoPrRed.^2).^0.5 ) * deltaZ;
+
+% Surface integral of dyadic product n*n (diagonal axisymmetric matrix):
+nn11 = pi * sum( rhoRed ./ (1 + rhoPrRed.^2).^0.5 ) * deltaZ; % = n22
+nn33 = 2*pi * sum( rhoRed .* rhoPrRed.^2 ./ (1 + rhoPrRed.^2).^0.5 ) * deltaZ;
+
+
+% --------------- Miscellaneous checks and plots --------------------------
+
+% Check whether surf = tr(integral n*n)
+if (surf - (2*nn11 + nn33) ) / surf > 0.0025 % strict equality beyond reach
+ disp(surf)
+ disp(2*nn11 + nn33)
+ startStr = ['r = ',num2str(rRatio),';theta = ',num2str(theta),';u* = ',num2str(uStar)];
+ error(['Liq-gas area # tr(integral over liq-gas surf. n*n) for ',startStr,';delta1 = ',num2str(delta1),';dZ = ',num2str(deltaZ)])
+end
+
+if plot==1 % && dist >=0 && vol >=0 %&& isfinite(vol)
+ plotProfile(rho,deltaZ,rRatio,delta1,delta2,theta,cstC,uStar)
+end
+
+if plot3D==1
+ plot3Dbridge(rho,deltaZ,rRatio)
+end
+
+% Check if the neck has been reached: (there is a risk because of divergence possibility)
+if uStar~= 0 % "y0" computation according to [10] Lian1993, (or (2.54) Soulie2005 - (11) Soulie2006)
+ rhoNeck = ( -1 + sqrt(1+2*uStar*cstC) ) / uStar;
+else
+ rhoNeck = cstC;
+end
+
+if (min(rho) - rhoNeck) / rhoNeck > 0.002
+ strU = num2str(uStar); strD1 = num2str(delta1); strDz = num2str(deltaZ);
+ if speak ~= 0
+ disp(['Profile with u*=',strU,', delta1=',strD1,' and deltaZ=',strDz,' has not reached its neck !'])
+ end
+ out = 0;
+else
+ out = 1;
+end
+
+end
+
+
+function angleRad = radians(angle)
+angleRad = angle * pi / 180;
+end
+
+function y = cstC_delta2(delta2,theta,uStar) % from e.g. (2.52) Soulie 2005
+global cstC
+y = sin(radians(delta2)) * sin(radians(delta2+theta)) + 1/2 * uStar * sin(radians(delta2))^2 - cstC;
+end
+
+
+function drho = drho(rho,prevRho,deltaZ,rho2d)
+% returns rhoPrime at i+1, see (4) Duriez2016
+
+%drho = sqrt( ( rho/(cstK - 1/2*uStar*rho^2) )^2 - 1 );% [11] Lian1993 always positiv, which is not true
+drho = (rho - prevRho) / deltaZ + 1/2*deltaZ * rho2d; % NB: the rho2d term has a real influence
+end
+
+function rho2d = rhoSecond(rho,rhoP,uStar)
+% see e.g. (5) Duriez2016
+rho2d = (1+rhoP^2) / rho + uStar*(1+rhoP^2)^1.5;
+end
+
+function plotProfile(rho,deltaZ,rRatio,delta1,delta2,theta,cstC,uStar)
+
+if uStar~= 0 % "y0" computation, see e.g. [10] Lian1993, (2.54) Soulie2005, (11) Soulie2006..
+ rhoNeck = ( -1 + sqrt(1+2*uStar*cstC) ) / uStar;
+else
+ rhoNeck = cstC;
+end
+lastZ = deltaZ*(length(rho) -1);
+vecZ = 0:deltaZ:lastZ;
+
+% Toroidal profile parameters: see [15]-[17] Lian1993, for r=1
+d1Rad = radians(delta1);d2Rad = radians(delta2); thRad = radians(theta);
+dist = lastZ - (1-cos(d1Rad)) - (1-cos(d2Rad)); % TODO: only for r=1
+rho1 = ( dist/2. + 1 - cos(d1Rad) )/cos(d2Rad+thRad);
+rho2 = sin(d2Rad) - (1-sin(d2Rad+thRad)) * ( dist/2 + 1-cos(d2Rad) ) / cos(d2Rad+thRad);
+
+figure();
+plot(vecZ,rho,'o');hold on;
+plot([0;lastZ],[1/rRatio * sin(radians(delta1));1/rRatio * sin(radians(delta1))],'r')
+plot([0;lastZ],[rhoNeck;rhoNeck],'m')
+plot([0;lastZ],[sin(radians(delta2));sin(radians(delta2))],'g')
+plot(vecZ,rho1 + abs(rho2) - ( rho1^2 - (vecZ-lastZ/2).^2 ).^(1/2),'b');
+legend('Profile','Left value','Neck value','Right value','Toroidal profile')
+xlabel('Interparticle axis, z')
+ylabel('\rho(z)')
+title(['\delta_1 = ',num2str(delta1)])
+if rRatio ~= 1
+ disp('Plotted toroidal profile is obviously wrong because rRatio # 1')
+end
+% disp([num2str(delta1),' profile plotted with ',num2str(length(rho)),' rho values'])
+
+end
+
+function plot3Dbridge(rho,deltaZ,rRatio)
+
+nPtsAlongAxis = 200 ; % grid-step along "rho" and z axis (the x,y plane of the 3D figure in fact)
+nRhoVal =length(rho);
+lastZ = deltaZ*(nRhoVal -1);
+x = -max(rho):max(rho)/nPtsAlongAxis:max(rho);
+
+y = 0:lastZ / (nPtsAlongAxis-1):lastZ; % will then reach exactly lastZ
+alt = zeros(length(y),length(x));
+for i = 1:length(y)
+ indRho = floor( (i-1)*(nRhoVal-1)/(nPtsAlongAxis) ) +1;
+ for j = 1:length(x)
+ if rho(indRho)^2 - x(j)^2 >=0
+ alt(i,j) = rho(indRho) * sqrt(rho(indRho)^2 - x(j)^2) / rho(indRho);
+% else
+% alt(i,j) = NaN; % not pretty... 2d step with matAltToPlot
+% necessary
+ end
+ end
+end
+
+figure();set(gca,'FontSize',14);
+surf(x,y,matAltToPlot(alt),'LineStyle','none') % filled surface, with color according to altitude
+% surf(x,y,matAltToPlot(alt),'FaceColor','none') % just the mesh
+ylabel('Interparticle axis, z')
+
+hold on
+
+delta1 = asin(rRatio*rho(1));
+yC1 = - cos(delta1)/rRatio;
+delta2 = asin(rho(nRhoVal));
+yC2 = lastZ + cos(delta2); % Center of right sphere:
+x = -max(rho):max(rho)/40:max(rho);
+y = x + lastZ/2;
+altSph1 = zeros(length(x),length(x));
+altSph2 = zeros(length(x),length(x));
+for i=1:length(x)
+ for j=1:length(x)
+ alt1 = 1/rRatio^2 - x(j)^2 - (y(i)-yC1)^2;
+ if alt1>=0
+ altSph1(i,j) = alt1^(1/2);
+ end
+ alt2 = 1 - x(j)^2 - (y(i)- yC2)^2;
+ if alt2>=0
+ altSph2(i,j) = alt2^(1/2);
+ end
+ end
+end
+
+surf(x,y,matAltToPlot(altSph1),'FaceColor','none')
+surf(x,y,matAltToPlot(altSph2),'FaceColor','none')
+caxis([0,max(rho)])
+grid off
+axis square
+end
+
+function retMat = matAltToPlot(matAlt)
+
+sizeY = size(matAlt,1);sizeX = size(matAlt,2);
+plotAlt = ones(sizeY,sizeX);
+for i = 2:sizeY-1
+ for j = 2:sizeX-1
+ if ~matAlt(i,j)
+ test = matAlt(i-1,j-1) || matAlt(i-1,j) || matAlt(i-1,j+1);
+ test = test || matAlt(i,j-1) || matAlt(i,j+1);
+ test = test || matAlt(i+1,j-1) || matAlt(i+1,j) || matAlt(i+1,j+1);
+ if ~test
+ plotAlt(i,j) = 0;
+ end
+ end
+ end
+end
+
+retMat = zeros(sizeY,sizeX);
+for i = 1:sizeY*sizeX
+ if ~plotAlt(i)
+ retMat(i) = NaN;
+ else
+ retMat(i) = matAlt(i);
+ end
+end
+
+
+end
+
=== added file 'examples/capillaryLaplaceYoung/writesCapFile.m'
--- examples/capillaryLaplaceYoung/writesCapFile.m 1970-01-01 00:00:00 +0000
+++ examples/capillaryLaplaceYoung/writesCapFile.m 2016-09-21 16:48:23 +0000
@@ -0,0 +1,254 @@
+% J. Duriez (jerome.duriez@xxxxxxxxxxx)
+
+function writesCapFile()
+% Builds capillary files for given radius ratio and contact angle
+% The historical template of YADE capillary files is kept: each line of
+% the capillary file is a bridge configuration. Configurations are
+% grouped by distance values, then (for any distance value) by capillary
+% pressure values.
+% The user is expected to modify l. 10 to 27 to taylor the generation to
+% his needs.
+
+% ---------------- Parameters to define by the user -----------------------
+
+% Laplace-Young resolution parameters:
+deltaZ = 2e-6; % use 2e-6, typically
+theta = 0; % contact angle, in degrees
+dDelta1 = 0.3; % delta1 (see solveLiqBridge.m) increment between two configurations
+
+% Capillary files construction parameters:
+nValUc = 350; % # (at most) of capillary pressure values to consider for
+% one given distance
+nValDist = 80; % # of distance values (linearly spaced between 0 and some
+% rupture distance computed in Preliminary A below) to consider
+% See also if needed l. 147 that defines the capillary file names
+
+% List of radius ratio to consider:
+rRatioVec = [1 1.1 1.25 1.5 1.75 2 3 4 5 10]; % historical rRatio values in YADE
+% As many capillary files as length(rRatioVec) will be generated
+
+suffixe = ''; % string to append at the end of generated files names
+
+% ---- End of parameters definitions, no further action is required ! -----
+
+
+
+
+% ---------------- Some data vector initializations -----------------------
+
+% maximum possible interparticle distances for given rRatio:
+maxD = -ones(length(rRatioVec),1);
+% Private (J. Duriez) remark: for theta=0; dz=1e-5 or 2e-6; delta1=15:0.2:90;
+% and rRatioVec = [1;1.1;1.25;1.5;1.75;2;3;4;5], we get here:
+% maxD = [0.3995;0.3808;0.3567;0.3241;0.2983;0.2772;0.2199;0.1852;0.1615];
+% Compares favourably with historical (L. Scholtes') data: dRupt_rMeSch.eps
+
+% initial guess (e.g. from L. Scholtes data) of maximum uc values for given rRatio:
+maxUc_r = [1 1493;1.1 1567;1.25 1781;1.5 2137;1.75 2469;2 2849;3 4274;4 5644;5 6987;10 14370];
+% these initial guess are not mandatory but may help save some time below
+
+
+% ---------- Generation of all capillary files really starts here ---------
+for i =1:length(rRatioVec)
+ rRatio = rRatioVec(i);
+
+% ##### Preliminary A: find maxD values #####
+% It is considered to be for zero suction (corresponds ~ to Scholtes' data)
+% And assumed to occur for delta1 greater than 15 deg (> 30 deg in Scholtes)
+ disp('Search for maximum distance starting')
+% fflush(stdout); % Octave command for flush display
+ d1here = 15:0.5:90-theta; % I got the exact same results with 15:1:90 or 15:0.2:90
+ data = solveLaplace_uc(theta,rRatio,0,d1here,deltaZ,0);
+ maxD(i) = max(data(:,1));
+ disp(['For r = ',num2str(rRatio),', maximum distance found: maxD = ',num2str(maxD(i))])
+ disp('');
+% fflush(stdout); % Octave command for flush display
+
+% #### Preliminary B: computing the set of d* values to consider, valDist ####
+
+ % -- Following lines are for geometric exponential sequence --
+ % valDist was a ~ 10-exponential sequence in L. Scholtes' data, but I
+ % (JD) think there was too much small values
+% valDist = -ones(nValDist,1);
+% valDist(1)=0;valDist(2) = 5e-5;% 1e-4 in Scholtes maybe too high (cf e.g. F* for r=3)
+% cst = 1 / (nValDist-2) * log(maxD(i) / valDist(2));
+% valDist(3:nValDist) = valDist(2) * exp( ((3:nValDist)'-2) * cst);
+
+ % -- Anyway, I (JD) think a linear sequence is just as good, or better --
+ incDist = maxD(i) / (nValDist - 1);
+ valDist = 0:incDist:maxD(i);
+ if length(valDist) ~= nValDist
+ error('We did not construct a consistent set of distances')
+ end
+
+% #### Preliminary C: find the greatest admissible suction for r,theta ####
+% It is considered to be the suction leading, for contacting particles, to
+% a mean filling angle <= 1 deg (was ~ 2 deg in L. Scholtes' data).
+
+ d1MaxUc = 0:dDelta1:min(45,90-theta);
+ disp('Starting search for maximum suction');%fflush(stdout);
+ if isempty(find(maxUc_r(:,1)==rRatio,1))
+ index = find(maxUc_r(:,1) < rRatio,1,'last');
+ ucTried = floor(maxUc_r(index,2));
+ else
+ ucTried = floor(maxUc_r(maxUc_r(:,1)==rRatio,2));
+ end
+ data = solveLaplace_uc(0,rRatio,ucTried,d1MaxUc,deltaZ,0);
+
+% Filling angles for contact situation are, if necessary, extrapolated:
+ deltaC = guessDeltaContact(data);
+ while deltaC > 1
+ prevMaxDelta1 = max(data(:,5));
+ d1MaxUc = 0:dDelta1:prevMaxDelta1;
+ ucTried = floor(ucTried*1.2); % floor because no need to take care of digits after comma here
+ data = solveLaplace_uc(0,rRatio,ucTried,d1MaxUc,deltaZ,0);
+ deltaC = guessDeltaContact(data);
+ end
+ maxUc = ucTried ;
+ disp(['For r = ',num2str(rRatio),', maximum suction found: final uc = ',num2str(maxUc)])
+ disp('');%fflush(stdout);
+% #### End of Preliminary C ####
+
+
+% #### Preliminary D: choosing list of uc* values, valUc #####
+% we define valUc as a geometric sequence --with 0 as first term, though--
+% so that log(uc*_{i+1}) - log(uc*_i) = cst, and reaching maxUc as the last
+% (i.e. nValUc-th) term
+ valUc = -ones(nValUc,1);
+ valUc(1)=0;valUc(2) = 0.002;% 0.002 as first positiv uc* value corresponds ~ to Scholtes
+ cst = 1 / (nValUc-2) * log(maxUc / valUc(2))/log(10);
+ valUc(3:nValUc) = valUc(2) * 10.^( ((3:nValUc)'-2) * cst);
+
+
+% #### Preliminaries completed: computing data for all uc values in valUc ! #####
+ bigData = -ones(nValDist*nValUc,8); % will gather all data
+ lineInBigData = 1;
+ disp('Starting solving for all uc values');tic
+% fflush(stdout);
+ output = zeros(10);iOut = 1; % these things are just to control some message display
+ for indUc = 1:nValUc
+ % message display:
+ if indUc >= iOut * nValUc / length(output) && output(iOut) == 0
+ disp([num2str(iOut / length(output)* 100),'% of uc values just completed'])
+% fflush(stdout);
+ toc;tic;
+ output(iOut) = 0;
+ iOut = iOut+1;
+ end
+ if indUc == 1
+ delta1Int = 0:dDelta1:90-theta;
+ end
+ smallData = solveLaplace_uc(theta,rRatio,valUc(indUc),delta1Int,deltaZ,0);
+ delta1Int = 0:dDelta1:max(smallData(:,5));
+ % Filling bigData:
+ nLinesSmallData = length(smallData(:,1));
+ bigData(lineInBigData:lineInBigData+nLinesSmallData-1,:) = smallData(:,[1:6,8:9]);
+ lineInBigData = lineInBigData + nLinesSmallData;
+ end
+ disp('All uc values computed !');toc
+% fflush(stdout);
+
+ % Cleaning and sorting the computed data:
+ bigData = bigData(bigData(:,1)>=0,:); % removal of possible extra data, due to preallocation
+ bigData = sortrows(bigData); % sorts according to increasing distance
+
+ nomCapFic = ['M(r=',num2str(rRatio),')',suffixe];
+ fic = fopen(nomCapFic,'w');
+ % first line of the file=rRatio value:
+ fprintf(fic,[num2str(rRatio),'\n',num2str(nValDist),'\n']);% for direct use in YADE
+ % alternatively, you may want to save some details about the
+ % construction. If yes, replace with the 2 following lines
+% fprintf(fic,[num2str(rRatio),' ; theta = ',num2str(theta),' ; dDelta1 = ',num2str(dDelta1)]);
+% fprintf(fic,[' ; deltaZ less or equal than ',num2str(deltaZ),'\n',num2str(nValDist),'\n']);
+ % in this case, you may need to manually suppress these extra
+ % characters before using the capillary file in YADE
+ fclose(fic);
+
+ % the final data have to conform valDist data. This is achieved through
+ % interpolation, below:
+ for indD =1:nValDist
+ dWanted = valDist(indD);
+ dataCapFile = interpolate(dWanted,valUc,bigData);
+ dlmwrite(nomCapFic,dataCapFile,'-append','delimiter','\t')
+ end
+
+
+
+end
+
+end
+
+function deltaC = guessDeltaContact(data)
+
+data = sortrows(data); % sorts the data according to increasing distance
+dist1 = data(1,1);dist2 = data(2,1); % the two lowest computed distances
+data1 = data(1,:); data2 = data(2,:); % the data for these two distances
+
+data_contact = data1 - ( data2 - data1 ) / (dist2 - dist1) * dist1;
+deltaC = mean(data_contact(5:6)); % mean delta1, delta2
+end
+
+function outData = interpolate(dist,valUc,bigData)
+% input bigData is a distance sorted bridges data table (matrix), obtained
+% directly from solving Laplace-Young equ. It hence shows no specific
+% distance values pattern.
+% Output outData is a data table for the given distance "dist" (and many uc
+% values) which is obtained distance-interpolating from bigData
+
+dataGreaterDist = bigData(bigData(:,1)>= dist,:);
+
+if isempty(dataGreaterDist)
+ % this should never happen: bigData includes all cases uc=0, and dist =
+ % maxD (the worst case) is one of them
+ outData = [];
+ out='\nUnexpected behavior obtained during interpolate function for dist =';
+ fprintf([out,num2str(dist),'. Please give a look..\n']);
+ return
+end
+maxUc_dist = max(dataGreaterDist(:,2));
+valUcOK = valUc(valUc <= maxUc_dist);
+nValUcOK = length(valUcOK);
+
+data_dist = -ones(nValUcOK,8); % all data for one given distance
+
+for i = 1:nValUcOK
+ currUc = valUcOK(i);
+ data_uc = bigData(bigData(:,2) == currUc,:);
+ if size(data_uc,1) == 1 % one single liq bridge data is useless
+ nValUcOK = i-1;
+ break
+ end
+ if dist < min(data_uc(:,1))
+ data1 = data_uc(1,:);
+ data2 = data_uc(2,:);
+ dist1 = data1(1) ; dist2 = data2(1);
+ data_dist(i,:) = data1 + ( data2 - data1 ) / (dist2 - dist1) * (dist - dist1);
+ else
+ data_dInf = data_uc(data_uc(:,1) <= dist,:); % not empty after previous if
+ data_dInf = data_dInf(size(data_dInf,1),:);
+ dInf = data_dInf(1);
+
+ data_dSup = data_uc(data_uc(:,1) >= dist,:); % not empty either
+ data_dSup = data_dSup(1,:);
+ dSup = data_dSup(1);
+
+ data_dist(i,:) = data_dInf + (data_dSup-data_dInf)/(dSup-dInf) * (dist-dInf);
+ end
+end
+
+data_dist = data_dist(data_dist(:,2) >= 0,:); % in case we faced the break above
+
+if size(data_dist,1) > nValUcOK
+ error('Problem here !')
+end
+if size(data_dist,1) == 0
+ % may happen for dist = maxD
+ out = '\nWARNING: Reaching the last distance value was not possible here. Before use in YADE, you need to';
+ out = [out,' modify by hand the 2d line (where nValDist appears) of this capillary file:\n'];
+ out = [out,'Just substract 1 to the nValDist value appearing therein\n'];
+ fprintf([out,'Sorry about that..\n\n']);
+end
+
+outData = [nValUcOK 0 0 0 0 0 0 0;data_dist];
+
+end
Follow ups
