sandbox/easystab/droplet_inclined_pipe.m

A droplet coming out a pipe under gravity effects

The aim of this code is to determine the shape of a droplet coming out a vertical or sligthly inclined pipe. We use Newton method in order to solve the nonlinear governing equations and Keller arc continuation method on the droplet volume in order to follow solution branches in the bifurcation plane.

The position of the interface is found solving the stress equilibrium at the interface. Similarly at the problem pendant_drop_volume.m, the stresses acting on the interface are given by the surface tension forces and the gravity which introduce an hydrostatic pressure inside the droplet. Thus we write the stress equilibrium as \displaystyle \rho g h - p = \sigma \kappa \, . where \sigma is the surface tension and \kappa is the interface curvature (here we consider only the curvature in the (x,y) plane).

Writing explicitly the curvature we obtain \displaystyle \rho g y - p = \sigma \frac{x^{''}y^{'}-y^{''}x^{'}}{(x^{'^2}+y^{'^2})^{3/2}} \, , where x y are the coordinates of the interface expressed in a parametric form. We decide to write the curvature in a parametric form in order to have a well posed problem also in the case where the interface is overturning (two y coordiantes of the interface corresponding to a single x coordinate).

Furthermore we impose conservation of volume \displaystyle \int_0^1 y x^' ds = V and we define the metrics of the parametric description of the interface \displaystyle x^{'^2}+y^{'^2}=l^2 where l is the length of the interface.

clear all; clf;

set(0,'defaultaxesfontsize',20,'defaultaxeslinewidth',.7,'defaultlinelinewidth',2,'defaultpatchlinewidth',.7);

    % parameters
    V = 2;                          % starting volume
    L = 1;                          % domain length
    n = 100;                        % number of grid points
    rho = 1;                        % density
    sigma = 1;                      % surface tension
    g = 1;                          % gravity
    theta = -0.0001;               % pipe inclination
    delta = 0.2;                    % continuation step

    % Differentiation and integration
    [D,DD,wt,t]=dif1D('cheb',0,1,n,3);
    Z=zeros(n,n); I=eye(n); 

Boundary conditions

Since we have a system containing a second and a first order differential equation we need three boundary conditions. Imposing the x and y coordinate for the starting and ending point of the interface we get four conditions. We use three of them as boundary conditions and we can still adopt the last one as fourth governing equation which close the system of four unknowns x, y, p, l.

    % Boundary conditions
    x1 = -L*cos(theta); xn = L*cos(theta);
    y1 = L*sin(theta); yn = -L*sin(theta);

Assuming a small volume and thus small interface deformations we can retrieve an intial guess by linearization which reads x = L(2t-1), y = \frac{p}{2\sigma}(L^2-x^2), p = \frac{3\sigma V}{2L^3} and l = \int_0^1 \sqrt{x'^2+y'^2}dt where 0<t<1 by definition.

    %initial guess
    p = 3*sigma*V/L^3/2;
    x = 2*t*L-L;
    y = p/2/sigma*(L^2-x.^2);
    l = sum(sqrt(diff(x).^2+diff(y).^2));

The solution vector contains the coordinates 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 in order to find different branches of the solution.

    sol = [x; y; p; l; V];
    dir=[zeros(2*n+2,1); 1];   % initial direction

In the following, 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.

    %initialization
    f = zeros(numel(sol),1);

    disp('Continuation loop')
    ind=1;quitcon=0;
    while ~quitcon

        solprev = sol;

Keller’s pseudo-arclength continuation

In order to follow different branches of the solution and turn around bifurcations we implement Keller’s pseudo-arclength continuation. For a more detailed explanation of the method look at Lecture Notes on Numerical Analysis of Nonlinear Equations, by Eusebius Doedel..

        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. See meniscus.m for more details. The implementation follows the one adopted in pendant_drop_volume.m where we find same governing equation and same type of boundary conditions. \displaystyle \rho g y - p - \sigma \frac{x^{''}y^{'}-y^{''}x^{'}}{(x^{'^2}+y^{'^2})^{3/2}} \, , the metric equation \displaystyle \quad x^{'^2}+y^{'^2}-l^2 \, , the constraint for the volume \displaystyle \int_0^1 y x^' ds - V \, , the equation for the remaning boundary condition \displaystyle x(0) - x_1 \, , and the additional equation for the continuation loop: \displaystyle \mathbf{t} \cdot ((V,H)_{new}^T-(V,H)_{old}^T) - \delta \, . The system 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
        hx = sol(1:n); hpx = D*hx; hppx = DD*hx; hy = sol(n+1:2*n); hpy = D*hy; hppy = DD*hy;
        a = hpx.^2+hpy.^2;  % this is comfortable because ita appears very often

        %current values of pressure, interface length and volume
        p = sol(end-2);
        l = sol(end-1);
        V = sol(end);

        % NONLINEAR FUNCTIONS
        %continuity of normal stress
        f(1:n) = p*a.^1.5-rho*g*hy.*a.^1.5-(hppx.*hpy-hpx.*hppy)*sigma;
        %metrics
        f(n+1:2*n) = hpx.^2+hpy.^2-l^2;
        %volume conservation CAREFUL WITH FLIPPED DROPLET, YOU NEED TO FLIP AROUND THE ORIGIN!!!
        f(2*n+1) = wt*(hy.*hpx)-V;
        %BC
        f(2*n+2) = hx(1)-x1;
        %continuation
        f(2*n+3) = dir'*(sol-solprev)-delta;

Computing the Jacobian

Since f is here a vector quantity, the Jacobian is a block matrix, whose blocks 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.

        % analytical jacobian
        A = [diag(3*p*a.^0.5.*hpx)*D-diag(sigma*hpy)*DD+diag(sigma*hppy)*D-diag(3*rho*g*hy.*hpx.*a.^0.5)*D ...
            diag(3*p*a.^0.5.*hpy)*D-diag(sigma*hppx)*D+diag(sigma*hpx)*DD-diag(a.^1.5*rho*g)-diag(3*rho*g*hy.*hpy.*a.^0.5)*D...
            a.^1.5 zeros(n,2); ... %from continuity of normal stress
            diag(2*hpx)*D diag(2*hpy)*D zeros(n,1) -2*l*ones(n,1) zeros(n,1); ...   %from metrics
            [-wt.*hpy' wt.*hpx' zeros(1,2) -1] + [-hy(1) zeros(1,n-2) hy(n) zeros(1,n+3)];...   %from volume conservation
            1 zeros(1,2*n+2); ...   %from BOUNDARY CONDITION that I use as equation
            dir']; %from continuation

The Jacobian has to be modified accordingly to the boundary conditions. Remember that one of them was already imposed in the construction of the Jacobian as an additional equation. The remaining three Dirichlet boundary conditions are: \displaystyle x(1)=xn \quad , \quad y(1) = y1 \quad , \quad y(1) = yn

        % Boundary conditions
        f([n n+1 2*n]) = [hx(n)-xn; hy(1)-y1; hy(n)-yn];
        A([n n+1 2*n],1:end) = [[I(n,:) zeros(1,n+2) 0];...
            [zeros(1,n) I(1,:) zeros(1,2) 0]; [zeros(1,n) I(n,:) zeros(1,2) 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);
        subplot(2,1,1)
        plot(hx,hy,'b-',[x1 xn],[y1 yn],'k',[x1-2*sin(theta) x1],[y1-2*abs(cos(theta)) y1],...
            'k',[xn-2*sin(theta) xn],[y1-2*abs(cos(theta)) yn],'k'); axis equal; drawnow;
        disp([num2str(count) '  ' num2str(res)]);
        if count>100||res>1e5; disp('no convergence'); delta = delta/2; break; end
        if res<1e-5; quit=1; disp('converged'); continue; end

The Newton step

Computation of the solution, it is obtained by solving a linear system.

        % Newton step
        sol = sol-A\f;

        count=count+1;
        end

        %this is to have the well defined value of the parameter that I
        %monitor
        xprime = D(1,:)*sol(1:n);
        yprime = D(1,:)*sol(n+1:2*n);
        if xprime>0&&yprime>0
            alpha = atan(yprime/xprime);
        elseif xprime>0&&yprime<0
            alpha = atan(yprime/xprime)+2*pi;
        elseif xprime<0&&yprime>0
            alpha = atan(yprime/xprime)+pi;
        elseif xprime<0&&yprime<0
            alpha = atan(yprime/xprime)+pi;
        end

The new direction for the continuation

The new direction for the continuation is found by:

        % New direction
        dir=A\[zeros(2*n+2,1);1]; % nouvelle direction
        dir=dir/norm(dir); % normalization
        ind = ind+1;

Bifurcation diagram and droplet shape

In our problem we choose the droplet volume V as homotopy parameter and we monitor the contact angle \alpha of the droplet when touching the pipe, this gives an indication on the bending of the interface. The forces acting on the system are surface tension and gravity, the first one trying to keep the droplet all together against the action of its own weigth. When the system is perfectly symmetric (not inclined pipe) the droplet is distribuited evenly respect to the pipe axis, in this case raising the droplet volume the interface will start to bend in order to better contrast the action of the weigth with surface tension forces, till reaching a point where a further bending is not enough anymore and a volume reduction is required leading to a bifurcation in our diagram. On the other hand when the system is asymmetric (even very slightly inclined pipe) the droplet mass is not distribuited evenly but it will be more on one side, from this physical observation we can figure out that for the same volume the required bending of the interface will be much higher that in the case before, thus the bifurcation will be reached for a smaller volume. The symmetry of this second phenomenon around the angle \theta will give rise to a symmetric bifurcation depending on this parameter (pitch-fork bifurcation). In the figure below we plot the contact angle \alpha of the droplet as a function of the volume as homotopy parameter, finding the bifurcations previously described. A saddle node bifuraction appears around V=5.2 if the pipe is not inclined (\theta=0) and a pitch-fork bifuraction around V=4.1 when the pipe is sligthly inclined (\theta=\pm 10^{-5}).

       subplot(2,1,2)
       hold on
       plot(V,alpha,'ob')
       hold off
       grid on
       title('bifurcation diagram')
       xlabel('volume')
       ylabel('\alpha')

       %break if it diverges
       if res>1e5
           break
       end
       
    end

The figure

Result of the computation.

Bifurcation diagram and the two possible droplet evolution

Exercices/Contributions

Play with the code trying different values of the inclination angle \theta, try to guess how it would look like the bifurcation diagram if we determine \alpha as a function of \theta. In order to study more in detail the influence of the angle \theta on this physical system one could use it as the homotopy parameter keeping the volume constant. Modify the code in this way and than try to use both the angle \theta and volume as continuation parameters. This will allow for a deeper understanding of the influnce of these parameters on the change in toplogy of the bifurcation diagram.