%{
# The shape of a pendant drop and the critical volume for dripping
The aim of this code is to compute the shape of a pendant drop situated
at the lower end of a vertical pipe, assumed plane. With the use of a
continuation procedure, the bifurcation diagram for increasing volume
of the droplet is also comptued.
Similarly as in the problem of the shape of the [meniscus.m](http://basilisk.fr/sandbox/easystab/meniscus.m),
the gravity induces a hydrostatic pressure distribution and the surface
tension imposes a pressure jump at the interface because of its
curvature. Furthermore, a pressure corresponds to the volume of the droplet.
The equation for the height of the interface $h$ reads:
$$
\sigma \frac{h^{''}}{(1+h^{'^2})^{3/2}}=\rho g h - p \, .
$$
For this problem so as for others where overturning phenomena are
present, it is convenient to consider a parametric description with the
interface being located at $x(s)$ and $y(s)$ in a cartesian coordinate
system with $dx^2+dy^2=l^2ds^2$ and $s\in[0,1]$, $l$ being the length of
the interface. The governing equation becomes then:
$$
\sigma \frac{x^{''}y^{'}-y^{''}x^{'}}{(x^{'^2}+y^{'^2})^{3/2}}=\rho g y - p \, .
$$
An additional equation is the conservation of the volume:
$$
\int_0^1 y x^' ds = V
$$
and because of the metric we also have:
$$
x^{'^2}+y^{'^2}=l^2
$$
%}
clear all; clf;
% parameters
L=1; % domain length
N=200; % number of points
rho=1; % density
sig=1; % surface tension
g=-1; % gravity
hi = 0.1; % thickness of uniform film
V=hi*L; % volume of uniform film
delta=0.2; % continuation length
% differentiation and integration
scale=-2/L;
[s,DM] = chebdif(N,2);
D=DM(:,:,1)*scale;
DD=DM(:,:,2)*scale^2;
s=(s-1)/scale;
Z=zeros(N,N); I=eye(N);
INT=([diff(s)',0]+[0,diff(s)'])/2;
%{
# Boundary conditions
Since we have a second order differential equation and two first order equations we
need four boundary conditions.
We impose Dirichlet boundary conditions, representing the fixed
positions of the droplet at the pipe wall.
Note that one condition will be used as additional equation in the system.
%}
% Boundary conditions
x1=0; xn=L;
y1=0; yn=0;
%{
As a initial guess for the Newton iteration it is sufficient to consider a uniform film of
thickness $h_i=V/L$ with $V$ being the volume and $L$ the width of the
pipe since we start with a small volume. The initial pressure is $p_i = V/L \rho g$.
%}
%initial guess
Pinit=V/L*rho*g;
initguess_x=s*L;
initguess_y=hi*ones(N,1);
linit=sum(sqrt(diff(initguess_x).^2+diff(initguess_y).^2));
%{
The solution vector contains the coordiantes of the interface, the
pressure, the length of the interface and the volume. The volume of the droplet is the
homotopy parameter, which is varied to obtain the bifurcation diagram.
Initially, the direction for the continuation is just along the volume.
%}
sol=[initguess_x;initguess_y;Pinit;linit;V];
dir=[zeros(N,1);zeros(N,1);0;0;1]; % initial direction
%{
In the next, the first outer loop is for the continuation over the volume $V$.
The inner loop is instead needed to find the steady-state solution
with a Newton method.
%}
disp('Continuation loop')
ind=1;quitcon=0;
while ~quitcon
solprev=sol;
%{
# Keller's pseudo-arclength continuation
To construct the bifurcation diagram we employ the property
that a steady state cannot appear or disappear.
There are different procedures for the continuation. The Keller's
pseudo-arclength continuation is used in this code. This method allows
the direction of the variation of the control
parameter to change sign and permits therefore to continuate the
bifurcation diagram after the fold.
Once the non-linear solution for a given parameter set is found, a guess
for the new one is made and a Netwon iteration is performed to find it.
The future value for the volume is not imposed by found again by the
Newton method and is indeed the last component of the residual function
$f$. The new solution (V,H) has to satisfy:
$$
\mathbf{t} \cdot ((V,H)_{new}^T-(V,H)_{old}^T) = \delta \, ,
$$
where $\mathbf{t}$ is the direction of the continuation and
$\delta$ the continuation length.
For more informations see the [Lecture Notes on Numerical Analysis of Nonlinear Equations, by Eusebius Doedel.](http://indy.cs.concordia.ca/auto/notes.pdf).
%}
sol=sol+dir*delta; % new prediction of solution
% Newton iterations
disp('Newton loop')
quit=0;count=0;
while ~quit
%{
# Newton loop
Now we are in the Newtoon loop to find the shape of the droplet with a specific volume (still unknow). See
[meniscus.m](http://basilisk.fr/sandbox/easystab/meniscus.m) for
more details.
The main difference is that here the function $f$ of which the root
has to be found is is a vector quantity whose components are the
governing equation
$$
\rho g y - p - \sigma
\frac{x^{''}y^{'}-y^{''}x^{'}}{(x^{'^2}+y^{'^2})^{3/2}} \, ,
$$
the metric equation
$$
\quad x^{'^2}+y^{'^2}-l^2 \, ,
$$
the constraint for the volume
$$
\int_0^1 y x^' ds - V \, ,
$$
the equation for one boundary
$$
x(0)=0 \, ,
$$
condition and the additional equation
for the continuation loop:
$$
\mathbf{t} \cdot ((V,H)_{new}^T-(V,H)_{old}^T) - \delta \, .
$$
It is therefore a function of the interface coordinates $x$, $y$,
the pressure $p$ and also of the length of the former, $l$, via the
curvilinear parameter.
*Notation*: the underscript `p` indicate a derivative w.r.t. `s`
%}
% the present solution and its derivatives
x=sol(1:N); y=sol(N+1:2*N); P=sol(2*N+1); l=sol(2*N+2); V=sol(end); xp=D*x; xpp=DD*x; yp=D*y; ypp=DD*y; a=xp.^2+yp.^2;
% nonlinear function
f=[-rho*g*y.*a.^1.5-sig.*(xpp.*yp-ypp.*xp)+P*a.^1.5; ...
a-l^2; ...
INT*(y.*xp)-V; ...
x(1)-x1; ...
dir'*(sol-solprev)-delta];
%{
# Computing the Jacobian
Since $f$ is here a vector quantity, the Jacobian is a bloc matrix,
whose blocs
are the discretized Jacobian matrices of the different components
w.r.t. the variables. They are found by perturbing all variables in the
considered equation with small parameters and neglecting the nonlinear
terms, in the same way as what described in [meniscus.m](http://basilisk.fr/sandbox/easystab/meniscus.m).
%}
% analytical jacobian
A=[ 3*P*diag(a.^0.5.*xp)*D-3*rho*g*diag(y.*xp.*a.^0.5)*D-sig*diag(yp)*DD+sig*diag(ypp)*D,...
3*P*diag(a.^0.5.*yp)*D-3*rho*g*diag(y.*yp.*a.^0.5)*D+sig*diag(xp)*DD-sig*diag(xpp)*D-rho*g*diag(a.^1.5),...
a.^1.5,zeros(N,1),zeros(N,1); ...
2*diag(xp)*D,2*diag(yp)*D,zeros(N,1),-2*l*ones(N,1),zeros(N,1); ...
y'.*I(N,:)-y'.*I(1,:)-INT.*yp',INT.*xp',0,0,-1;...
I(1,:),zeros(1,N),0,0,0; ...
dir'];
%{
The Jacobian has to be modified accordingly to the boundary conditions as well.
Remember that one of them was already imposed as a
constraint in the construction of the Jacobian as an additional equation.
The remaining three Dirichlet boundary conditions are:
$$
x(1)=x_n \quad , \quad y(0) = y_1 \quad , \quad y(1) = y_n
$$
%}
% Boundary conditions
loc = [1 N N+1];
f(loc)=[x(N)-xn; y(1)-y1; y(N)-yn];
A(loc,:)=[I(N,:),zeros(1,N),0,0,0; zeros(1,N),I(1,:),0,0,0; zeros(1,N),I(N,:),0,0,0];
%{
If the residual is smaller than the tolerance, the shape of the droplet
for this volume is found. If this is not the case, an additional
Newton iteration has to be performed.
%}
% convergence test
res=norm(f);
%plot(x,y,'b-',initguess_x,initguess_y,'r--'); drawnow;
disp([num2str(count) ' ' num2str(res)]);
if count>50|res>1e5; disp('no convergence'); endLoop =1; quitcon=1; break; end
if res<1e-5; quit=1; disp('converged'); continue; end
%{
# The Newton step
The core of the Netwon iteration is the computation of the new guess
for the solution.
%}
% Newton step
sol=sol-A\f;
count=count+1;
end
%{
# The new direction for the continuation
The new direction for the continuation is found by:
%}
% New direction
dir=A\[zeros(N,1);zeros(N,1);0;0;1]; % new direction
dir=dir/max(abs(dir)); % normalization
ind=ind+1;
%{
# Bifurcation diagram and droplet shape
The shape of the pendant drop for a given volume is finally found and its
height can be plotted as a function of its volume to construct the bifurcation diagram.
A saddle node bifurcation is clearly visible from the diagram. After
the fold, the branch is unstable and the volume has to be reduced to
find a steady-state solution. Note that the stability of the solution cannot be
determined with this method. Here the stability character of the branches is
guessed from physical arguments.
%}
x=sol(1:N);
y=sol(N+1:2*N);
V=sol(end);
H = max(y);
figure(1)
plot(x,-y,'b-',[0 L],[0 0],'-k'); axis equal; %axis([-L/2 3/2*L 0 4*L]);
figure(2)
plot(V,H,'bo')
hold on
grid on
xlabel('V')
ylabel('H')
drawnow;
end
%{
# The figure
Here we see the result of the computation.
![Bifurcation diagram and the different shapes of droplet](/pendant_drop_volume.png)
%}
%{
# Exercices/Contributions
Modify the previous code to consider the problem of a pendant droplet connected
to a uniform film (zero contact angle) or to a horizontal substrate ( given contact angle depending on materials).
This is for example the situation of the dripping from
the ceiling.
Start by implementing Neumann boundary conditions to
take into account the inclination of the interface at its extremities.
%}