Vous avez reçu un message "Your GitLab account has been locked ..." ? Pas d'inquiétude : lisez cet article https://docs.gricad-pages.univ-grenoble-alpes.fr/help/unlock/

Commit e23a181f authored by jduriez's avatar jduriez
Browse files

Capillary scripts commit...

Capillary scripts commit (http://www.mail-archive.com/yade-dev@lists.launchpad.net/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
parent 8408148b
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.
# 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)
% J. Duriez (jerome.duriez@ucalgary.ca)
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
% J. Duriez (jerome.duriez@ucalgary.ca)
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
% J. Duriez (jerome.duriez@ucalgary.ca)
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);